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/UKMO/NEMO_4.0-TRUNK_r14960_HPG/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0-TRUNK_r14960_HPG/src/OCE/DYN/dynhpg.F90 @ 15726

Last change on this file since 15726 was 15726, checked in by dbruciaferri, 2 years ago
adding V1 of DJR and FFR code - dynhpg_djr_ffr_smooth_v1.F90
Property svn:keywords set to `Id`
File size: 171.1 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 constrained cubic splines (ccs))
31	!! hpg_prj : s-coordinate (Pressure Jacobian with Cubic polynomial)
32	!! hpg_djr : s-coordinate (Density Jacobian with ccs subtracting a reference)
33	!! hpg_ffr : s-coordinate (Forces on faces subtracting a reference profile)
34	!!----------------------------------------------------------------------
35	USE oce ! ocean dynamics and tracers
36	USE isf_oce , ONLY : risfload ! ice shelf (risfload variable)
37	USE isfload , ONLY : isf_load ! ice shelf (isf_load routine )
38	USE sbc_oce ! surface variable (only for the flag with ice shelf)
39	USE dom_oce ! ocean space and time domain
40	USE wet_dry ! wetting and drying
41	USE phycst ! physical constants
42	USE trd_oce ! trends: ocean variables
43	USE trddyn ! trend manager: dynamics
44	USE zpshde ! partial step: hor. derivative (zps_hde routine)
45	!
46	USE in_out_manager ! I/O manager
47	USE prtctl ! Print control
48	USE lbclnk ! lateral boundary condition
49	USE lib_mpp ! MPP library
50	USE eosbn2 ! compute density
51	USE timing ! Timing
52	USE iom
53
54	IMPLICIT NONE
55	PRIVATE
56
57	PUBLIC dyn_hpg ! routine called by step module
58	PUBLIC dyn_hpg_init ! routine called by opa module
59
60	! !!* Namelist namdyn_hpg : hydrostatic pressure gradient
61	LOGICAL, PUBLIC :: ln_hpg_zco !: z-coordinate - full steps
62	LOGICAL, PUBLIC :: ln_hpg_zps !: z-coordinate - partial steps (interpolation)
63	LOGICAL, PUBLIC :: ln_hpg_sco !: s-coordinate (standard jacobian formulation)
64	LOGICAL, PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial)
65	LOGICAL, PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme)
66	LOGICAL, PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf
67	LOGICAL, PUBLIC :: ln_hpg_djr !: s-coordinate (density Jacobian with cubic polynomial, subtracting a local reference profile)
68	LOGICAL, PUBLIC :: ln_hpg_ffr !: s-coordinate (forces on faces with subtraction of a local reference profile)
69
70	! !! Flag to control the type of hydrostatic pressure gradient
71	INTEGER, PARAMETER :: np_ERROR =-10 ! error in specification of lateral diffusion
72	INTEGER, PARAMETER :: np_zco = 0 ! z-coordinate - full steps
73	INTEGER, PARAMETER :: np_zps = 1 ! z-coordinate - partial steps (interpolation)
74	INTEGER, PARAMETER :: np_sco = 2 ! s-coordinate (standard jacobian formulation)
75	INTEGER, PARAMETER :: np_djc = 3 ! s-coordinate (Density Jacobian with Cubic polynomial)
76	INTEGER, PARAMETER :: np_prj = 4 ! s-coordinate (Pressure Jacobian scheme)
77	INTEGER, PARAMETER :: np_isf = 5 ! s-coordinate similar to sco modify for isf
78	INTEGER, PARAMETER :: np_djr = 6 ! s-coordinate (density Jacobian with cubic polynomial, subtracting a local reference profile)
79	INTEGER, PARAMETER :: np_ffr = 7 ! s-coordinate (forces on faces with subtraction of a local reference profile)
80
81	!
82	INTEGER, PUBLIC :: nhpg !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM)
83	!
84
85	LOGICAL :: ln_hpg_bcvN_rhd_hor, ln_hpg_bcvN_rhd_srf ! flags to specify constrained cubic spline (ccs) bdy conditions for rhd & z
86	LOGICAL :: ln_hpg_bcvN_rhd_bot, ln_hpg_bcvN_z_hor ! True implies von Neumann bcs; False implies linear extrapolation
87
88	REAL(wp) :: aco_bc_rhd_hor, bco_bc_rhd_hor, aco_bc_rhd_srf ! coefficients for hpg_djc & hpg_djr boundary conditions
89	REAL(wp) :: bco_bc_rhd_srf, aco_bc_rhd_bot, bco_bc_rhd_bot ! "
90	REAL(wp) :: aco_bc_z_hor, bco_bc_z_hor, aco_bc_z_srf ! "
91	REAL(wp) :: bco_bc_z_srf, aco_bc_z_bot, bco_bc_z_bot ! "
92
93	LOGICAL :: ln_hpg_djr_ref_ccs ! T => constrained cubic, F => simple cubic used for interpolation of reference to target profiles
94	LOGICAL :: ln_hpg_ffr_ref ! T => reference profile is subtracted; F => no reference profile subtracted
95	LOGICAL :: ln_hpg_ffr_ref_ccs ! T => use constrained cubic spline to vertically interpolate reference to each profile
96	LOGICAL :: ln_hpg_ffr_hor_ccs ! T => use constrained cubic spline to interpolate z_rhd_pmr along upper faces of cells
97	LOGICAL :: ln_hpg_ffr_hor_cub ! T => use simple cubic polynomial to interpolate z_rhd_pmr along upper faces of cells
98	LOGICAL :: ln_hpg_ffr_vrt_quad ! T => use quadratic fit to z_rhd_pmr in vertical interpoln & integration; else standard 2nd order scheme
99
100	LOGICAL :: ln_dbg_hpg ! T => debug outputs generated
101	LOGICAL :: ln_dbg_ij, ln_dbg_ik, ln_dbg_jk
102	INTEGER :: ki_dbg_min, ki_dbg_max, ki_dbg_cen
103	INTEGER :: kj_dbg_min, kj_dbg_max, kj_dbg_cen
104	INTEGER :: kk_dbg_min, kk_dbg_max, kk_dbg_cen
105
106
107	!! * Substitutions
108	# include "do_loop_substitute.h90"
109	# include "domzgr_substitute.h90"
110
111	!!----------------------------------------------------------------------
112	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
113	!! $Id$
114	!! Software governed by the CeCILL license (see ./LICENSE)
115	!!----------------------------------------------------------------------
116	CONTAINS
117
118	SUBROUTINE dyn_hpg( kt, Kmm, puu, pvv, Krhs )
119	!!---------------------------------------------------------------------
120	!! * ROUTINE dyn_hpg *
121	!!
122	!! ** Method : Call the hydrostatic pressure gradient routine
123	!! using the scheme defined in the namelist
124	!!
125	!! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
126	!! - send trends to trd_dyn for futher diagnostics (l_trddyn=T)
127	!!----------------------------------------------------------------------
128	INTEGER , INTENT( in ) :: kt ! ocean time-step index
129	INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
130	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
131	!
132	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv
133	!!----------------------------------------------------------------------
134	!
135	IF( ln_timing ) CALL timing_start('dyn_hpg')
136	!
137	IF( l_trddyn ) THEN ! Temporary saving of puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends (l_trddyn)
138	ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
139	ztrdu(:,:,:) = puu(:,:,:,Krhs)
140	ztrdv(:,:,:) = pvv(:,:,:,Krhs)
141	ENDIF
142	!
143	SELECT CASE ( nhpg ) ! Hydrostatic pressure gradient computation
144	CASE ( np_zco ) ; CALL hpg_zco ( kt, Kmm, puu, pvv, Krhs ) ! z-coordinate
145	CASE ( np_zps ) ; CALL hpg_zps ( kt, Kmm, puu, pvv, Krhs ) ! z-coordinate plus partial steps (interpolation)
146	CASE ( np_sco ) ; CALL hpg_sco ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (standard jacobian formulation)
147	CASE ( np_djc ) ; CALL hpg_djc ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (Density Jacobian with Cubic polynomial)
148	CASE ( np_prj ) ; CALL hpg_prj ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (Pressure Jacobian scheme)
149	CASE ( np_isf ) ; CALL hpg_isf ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate similar to sco modify for ice shelf
150	CASE ( np_djr ) ; CALL hpg_djr ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (Density Jacobian with Cubic polynomial subtracting a local reference profile)
151	CASE ( np_ffr ) ; CALL hpg_ffr ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (forces on faces with subtraction of a local reference profile)
152	END SELECT
153	!
154	IF( l_trddyn ) THEN ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
155	ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
156	ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
157	CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt, Kmm )
158	DEALLOCATE( ztrdu , ztrdv )
159	ENDIF
160	!
161	IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' hpg - Ua: ', mask1=umask, &
162	& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
163	!
164	IF( ln_timing ) CALL timing_stop('dyn_hpg')
165	!
166	END SUBROUTINE dyn_hpg
167
168
169	SUBROUTINE dyn_hpg_init( Kmm )
170	!!----------------------------------------------------------------------
171	!! * ROUTINE dyn_hpg_init *
172	!!
173	!! ** Purpose : initializations for the hydrostatic pressure gradient
174	!! computation and consistency control
175	!!
176	!! ** Action : Read the namelist namdyn_hpg and check the consistency
177	!! with the type of vertical coordinate used (zco, zps, sco)
178	!!----------------------------------------------------------------------
179	INTEGER, INTENT( in ) :: Kmm ! ocean time level index
180	!
181	INTEGER :: ioptio = 0 ! temporary integer
182	INTEGER :: ios ! Local integer output status for namelist read
183	!!
184	INTEGER :: ji, jj, jk, ikt ! dummy loop indices ISF
185	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zts_top, zrhd ! hypothesys on isf density
186	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zrhdtop_isf ! density at bottom of ISF
187	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ziceload ! density at bottom of ISF
188	!!
189
190	NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_isf, &
191	& ln_hpg_djc, ln_hpg_prj, ln_hpg_djr, ln_hpg_ffr, &
192	& ln_hpg_bcvN_rhd_hor, ln_hpg_bcvN_rhd_srf, &
193	& ln_hpg_bcvN_rhd_bot, ln_hpg_bcvN_z_hor, &
194	& ln_hpg_djr_ref_ccs, ln_hpg_ffr_vrt_quad, &
195	& ln_hpg_ffr_ref, ln_hpg_ffr_ref_ccs, &
196	& ln_hpg_ffr_hor_ccs, ln_hpg_ffr_hor_cub, &
197	& ln_dbg_hpg, ln_dbg_ij, ln_dbg_ik, ln_dbg_jk, &
198	& ki_dbg_min, ki_dbg_max, ki_dbg_cen, &
199	& kj_dbg_min, kj_dbg_max, kj_dbg_cen, &
200	& kk_dbg_min, kk_dbg_max, kk_dbg_cen
201
202	!!----------------------------------------------------------------------
203	!
204	READ ( numnam_ref, namdyn_hpg, IOSTAT = ios, ERR = 901)
205	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist' )
206	!
207	READ ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 )
208	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist' )
209	IF(lwm) WRITE ( numond, namdyn_hpg )
210	!
211	IF(lwp) THEN ! Control print
212	WRITE(numout,*)
213	WRITE(numout,*) 'dyn_hpg_init : hydrostatic pressure gradient initialisation'
214	WRITE(numout,*) '~~~~~~~~~~~~'
215	WRITE(numout,*) ' Namelist namdyn_hpg : choice of hpg scheme'
216	WRITE(numout,*) ' z-coord. - full steps ln_hpg_zco = ', ln_hpg_zco
217	WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps
218	WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco
219	WRITE(numout,*) ' s-coord. (standard jacobian formulation) for isf ln_hpg_isf = ', ln_hpg_isf
220	WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc
221	WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj
222	WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic minus reference) ln_hpg_djr = ', ln_hpg_djr
223	WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic minus reference) ln_hpg_ffr = ', ln_hpg_djr
224	WRITE(numout,*) ' s-coord. (forces on faces minus local reference) ln_hpg_ffr = ', ln_hpg_ffr
225	ENDIF
226	!
227	IF( .NOT.ln_linssh .AND. (ln_hpg_zco.OR.ln_hpg_zps) ) &
228	& CALL ctl_stop( 'dyn_hpg_init : non-linear free surface incompatible with hpg_zco or hpg_zps' )
229	!
230	IF( (.NOT.ln_hpg_isf .AND. ln_isfcav) .OR. (ln_hpg_isf .AND. .NOT.ln_isfcav) ) &
231	& CALL ctl_stop( 'dyn_hpg_init : ln_hpg_isf=T requires ln_isfcav=T and vice versa' )
232
233	!
234	#if defined key_qco
235	IF( ln_hpg_isf ) THEN
236	CALL ctl_stop( 'dyn_hpg_init : key_qco and ln_hpg_isf not yet compatible' )
237	ENDIF
238	#endif
239
240	IF( ln_hpg_djr .AND. nn_hls < 2) THEN
241	CALL ctl_stop( 'dyn_hpg_init : nn_hls < 2 and ln_hpg_djr not yet compatible' )
242	ENDIF
243
244	!
245	! ! Set nhpg from ln_hpg_... flags & consistency check
246	nhpg = np_ERROR
247	ioptio = 0
248	IF( ln_hpg_zco ) THEN ; nhpg = np_zco ; ioptio = ioptio +1 ; ENDIF
249	IF( ln_hpg_zps ) THEN ; nhpg = np_zps ; ioptio = ioptio +1 ; ENDIF
250	IF( ln_hpg_sco ) THEN ; nhpg = np_sco ; ioptio = ioptio +1 ; ENDIF
251	IF( ln_hpg_djc ) THEN ; nhpg = np_djc ; ioptio = ioptio +1 ; ENDIF
252	IF( ln_hpg_prj ) THEN ; nhpg = np_prj ; ioptio = ioptio +1 ; ENDIF
253	IF( ln_hpg_isf ) THEN ; nhpg = np_isf ; ioptio = ioptio +1 ; ENDIF
254	IF( ln_hpg_djr ) THEN ; nhpg = np_djr ; ioptio = ioptio +1 ; ENDIF
255	IF( ln_hpg_ffr ) THEN ; nhpg = np_ffr ; ioptio = ioptio +1 ; ENDIF
256	!
257	IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
258	!
259	IF(lwp) THEN
260	WRITE(numout,*)
261	SELECT CASE( nhpg )
262	CASE( np_zco ) ; WRITE(numout,*) ' ==>>> z-coord. - full steps '
263	CASE( np_zps ) ; WRITE(numout,*) ' ==>>> z-coord. - partial steps (interpolation)'
264	CASE( np_sco ) ; WRITE(numout,*) ' ==>>> s-coord. (standard jacobian formulation)'
265	CASE( np_djc ) ; WRITE(numout,*) ' ==>>> s-coord. (Density Jacobian: Cubic polynomial)'
266	CASE( np_prj ) ; WRITE(numout,*) ' ==>>> s-coord. (Pressure Jacobian: Cubic polynomial)'
267	CASE( np_isf ) ; WRITE(numout,*) ' ==>>> s-coord. (standard jacobian formulation) for isf'
268	CASE( np_djr ) ; WRITE(numout,*) ' ==>>> s-coord. (Density Jacobian: Cubic polynomial subtracting a local reference profile)'
269	CASE( np_ffr ) ; WRITE(numout,*) ' ==>>> s-coord. (forces on faces subtracting a local reference profile)'
270
271	END SELECT
272	WRITE(numout,*)
273	ENDIF
274	!
275	IF ( ln_hpg_djc .OR. ln_hpg_djr .OR. ln_hpg_ffr ) THEN
276	IF (ln_hpg_bcvN_rhd_hor) THEN ! Von Neumann boundary condition
277	IF(lwp) WRITE(numout,*) ' ccs rhd horizontal bc: von Neumann '
278	aco_bc_rhd_hor = 6.0_wp/5.0_wp
279	bco_bc_rhd_hor = 7.0_wp/15.0_wp
280	ELSE ! Linear extrapolation
281	IF(lwp) WRITE(numout,*) ' ccs rhd horizontal bc: linear extrapolation'
282	aco_bc_rhd_hor = 3.0_wp/2.0_wp
283	bco_bc_rhd_hor = 1.0_wp/2.0_wp
284	END IF
285	IF (ln_hpg_bcvN_rhd_srf) THEN ! Von Neumann boundary condition
286	IF(lwp) WRITE(numout,*) ' ccs rhd surface bc: von Neumann '
287	aco_bc_rhd_srf = 6.0_wp/5.0_wp
288	bco_bc_rhd_srf = 7.0_wp/15.0_wp
289	ELSE ! Linear extrapolation
290	IF(lwp) WRITE(numout,*) ' ccs rhd surface bc: linear extrapolation'
291	aco_bc_rhd_srf = 3.0_wp/2.0_wp
292	bco_bc_rhd_srf = 1.0_wp/2.0_wp
293	END IF
294	IF (ln_hpg_bcvN_rhd_bot) THEN ! Von Neumann boundary condition
295	IF(lwp) WRITE(numout,*) ' ccs rhd bottom bc: von Neumann '
296	aco_bc_rhd_bot = 6.0_wp/5.0_wp
297	bco_bc_rhd_bot = 7.0_wp/15.0_wp
298	ELSE ! Linear extrapolation
299	IF(lwp) WRITE(numout,*) ' ccs rhd bottom bc: linear extrapolation'
300	aco_bc_rhd_bot = 3.0_wp/2.0_wp
301	bco_bc_rhd_bot = 1.0_wp/2.0_wp
302	END IF
303	IF (ln_hpg_bcvN_z_hor) THEN ! Von Neumann boundary condition
304	IF(lwp) WRITE(numout,*) ' ccs z horizontal bc: von Neumann '
305	aco_bc_z_hor = 6.0_wp/5.0_wp
306	bco_bc_z_hor = 7.0_wp/15.0_wp
307	ELSE ! Linear extrapolation
308	IF(lwp) WRITE(numout,*) ' ccs z horizontal bc: linear extrapolation'
309	aco_bc_z_hor = 3.0_wp/2.0_wp
310	bco_bc_z_hor = 1.0_wp/2.0_wp
311	END IF
312
313	! linear extrapolation used for z in the vertical (surface and bottom)
314	aco_bc_z_srf = 3.0_wp/2.0_wp
315	bco_bc_z_srf = 1.0_wp/2.0_wp
316	aco_bc_z_bot = 3.0_wp/2.0_wp
317	bco_bc_z_bot = 1.0_wp/2.0_wp
318
319	END IF
320
321	IF ( lwp .AND. ln_hpg_djr ) THEN
322	IF ( ln_hpg_djr_ref_ccs ) THEN
323	WRITE(numout,*) ' interpolation of reference profile by constrained cubic spline'
324	ELSE
325	WRITE(numout,*) ' interpolation of reference profile by simple cubic spline with off-centring at boundaries'
326	END IF
327	END IF
328
329	IF ( lwp .AND. ln_hpg_ffr ) THEN
330	IF ( ln_hpg_ffr_ref_ccs ) THEN
331	IF ( ln_hpg_ffr_ref_ccs ) THEN
332	WRITE(numout,*) ' interpolation of reference profile in vertical by constrained cubic spline '
333	ELSE
334	WRITE(numout,*) ' interpolation of reference profile in vertical by pure cubic polynomial fit '
335	END IF
336	ELSE
337	WRITE(numout,*) ' no reference profile subtracted '
338	END IF
339	IF ( ln_hpg_ffr_hor_ccs ) THEN
340	WRITE(numout,*) ' interpolation of z_rhd_pmr profile in horizontal by constrained cubic spline '
341	ELSE if ( ln_hpg_ffr_hor_cub ) THEN
342	WRITE(numout,*) ' interpolation of z_rhd_pmr profile in horizontal by simple cubic polynomial '
343	ELSE
344	WRITE(numout,*) ' simple linear interpolation in horizontal of z_rhd_pmr profile'
345	END IF
346	IF ( ln_hpg_ffr_vrt_quad ) THEN
347	WRITE(numout,*) ' quadratic fit to z_rhd_pmr profile in vertical for interpolation and integration '
348	ELSE
349	WRITE(numout,*) ' simple second order accurate vertical integration of z_rhd_pmr profile in vertical'
350	END IF
351
352	END IF
353
354	!
355	IF ( ln_dbg_hpg .AND. lwp ) THEN
356	WRITE(numout,*) ' dyn_hpg diagnostic settings'
357	WRITE(numout,*) ' ki_dbg_min = ', ki_dbg_min, '; ki_dbg_max = ', ki_dbg_max
358	WRITE(numout,*) ' kj_dbg_min = ', kj_dbg_min, '; kj_dbg_max = ', kj_dbg_max
359	WRITE(numout,*) ' kk_dbg_min = ', kk_dbg_min, '; kk_dbg_max = ', kk_dbg_max
360	WRITE(numout,*) ' ki_dbg_cen = ', ki_dbg_cen, '; kj_dbg_cen = ', kj_dbg_cen
361	WRITE(numout,*) ' kk_dbg_cen = ', kk_dbg_cen
362	END IF
363
364	END SUBROUTINE dyn_hpg_init
365
366
367	SUBROUTINE hpg_zco( kt, Kmm, puu, pvv, Krhs )
368	!!---------------------------------------------------------------------
369	!! * ROUTINE hpg_zco *
370	!!
371	!! ** Method : z-coordinate case, levels are horizontal surfaces.
372	!! The now hydrostatic pressure gradient at a given level, jk,
373	!! is computed by taking the vertical integral of the in-situ
374	!! density gradient along the model level from the suface to that
375	!! level: zhpi = grav .....
376	!! zhpj = grav .....
377	!! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
378	!! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi
379	!! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj
380	!!
381	!! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
382	!!----------------------------------------------------------------------
383	INTEGER , INTENT( in ) :: kt ! ocean time-step index
384	INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
385	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
386	!
387	INTEGER :: ji, jj, jk ! dummy loop indices
388	REAL(wp) :: zcoef0, zcoef1 ! temporary scalars
389	REAL(wp), DIMENSION(A2D(nn_hls)) :: zhpi, zhpj
390	!!----------------------------------------------------------------------
391	!
392	IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
393	IF( kt == nit000 ) THEN
394	IF(lwp) WRITE(numout,*)
395	IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend'
396	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case '
397	ENDIF
398	ENDIF
399	!
400	zcoef0 = - grav * 0.5_wp ! Local constant initialization
401	!
402	DO_2D( 0, 0, 0, 0 ) ! Surface value
403	zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm)
404	! ! hydrostatic pressure gradient
405	zhpi(ji,jj) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
406	zhpj(ji,jj) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
407	! ! add to the general momentum trend
408	puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj)
409	pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj)
410	END_2D
411	!
412	DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! interior value (2=<jk=<jpkm1)
413	zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm)
414	! ! hydrostatic pressure gradient
415	zhpi(ji,jj) = zhpi(ji,jj) + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) &
416	& - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj)
417
418	zhpj(ji,jj) = zhpj(ji,jj) + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) &
419	& - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj)
420	! ! add to the general momentum trend
421	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj)
422	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj)
423	END_3D
424	!
425	END SUBROUTINE hpg_zco
426
427
428	SUBROUTINE hpg_zps( kt, Kmm, puu, pvv, Krhs )
429	!!---------------------------------------------------------------------
430	!! * ROUTINE hpg_zps *
431	!!
432	!! ** Method : z-coordinate plus partial steps case. blahblah...
433	!!
434	!! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
435	!!----------------------------------------------------------------------
436	INTEGER , INTENT( in ) :: kt ! ocean time-step index
437	INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
438	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
439	!!
440	INTEGER :: ji, jj, jk ! dummy loop indices
441	INTEGER :: iku, ikv ! temporary integers
442	REAL(wp) :: zcoef0, zcoef1, zcoef2, zcoef3 ! temporary scalars
443	REAL(wp), DIMENSION(A2D(nn_hls),jpk ) :: zhpi, zhpj
444	REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: zgtsu, zgtsv
445	REAL(wp), DIMENSION(A2D(nn_hls) ) :: zgru, zgrv
446	!!----------------------------------------------------------------------
447	!
448	IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
449	IF( kt == nit000 ) THEN
450	IF(lwp) WRITE(numout,*)
451	IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend'
452	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate with partial steps - vector optimization'
453	ENDIF
454	ENDIF
455
456	! Partial steps: Compute NOW horizontal gradient of t, s, rd at the last ocean level
457	CALL zps_hde( kt, Kmm, jpts, ts(:,:,:,:,Kmm), zgtsu, zgtsv, rhd, zgru , zgrv )
458
459	! Local constant initialization
460	zcoef0 = - grav * 0.5_wp
461
462	! Surface value (also valid in partial step case)
463	DO_2D( 0, 0, 0, 0 )
464	zcoef1 = zcoef0 * e3w(ji,jj,1,Kmm)
465	! hydrostatic pressure gradient
466	zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj)
467	zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj)
468	! add to the general momentum trend
469	puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1)
470	pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1)
471	END_2D
472
473	! interior value (2=<jk=<jpkm1)
474	DO_3D( 0, 0, 0, 0, 2, jpkm1 )
475	zcoef1 = zcoef0 * e3w(ji,jj,jk,Kmm)
476	! hydrostatic pressure gradient
477	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &
478	& + zcoef1 * ( ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) &
479	& - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj)
480
481	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &
482	& + zcoef1 * ( ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) &
483	& - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj)
484	! add to the general momentum trend
485	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk)
486	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk)
487	END_3D
488
489	! partial steps correction at the last level (use zgru & zgrv computed in zpshde.F90)
490	DO_2D( 0, 0, 0, 0 )
491	iku = mbku(ji,jj)
492	ikv = mbkv(ji,jj)
493	zcoef2 = zcoef0 * MIN( e3w(ji,jj,iku,Kmm), e3w(ji+1,jj ,iku,Kmm) )
494	zcoef3 = zcoef0 * MIN( e3w(ji,jj,ikv,Kmm), e3w(ji ,jj+1,ikv,Kmm) )
495	IF( iku > 1 ) THEN ! on i-direction (level 2 or more)
496	puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) - zhpi(ji,jj,iku) ! subtract old value
497	zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one
498	& + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + zgru(ji,jj) ) * r1_e1u(ji,jj)
499	puu (ji,jj,iku,Krhs) = puu(ji,jj,iku,Krhs) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend
500	ENDIF
501	IF( ikv > 1 ) THEN ! on j-direction (level 2 or more)
502	pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) - zhpj(ji,jj,ikv) ! subtract old value
503	zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one
504	& + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + zgrv(ji,jj) ) * r1_e2v(ji,jj)
505	pvv (ji,jj,ikv,Krhs) = pvv(ji,jj,ikv,Krhs) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend
506	ENDIF
507	END_2D
508	!
509	END SUBROUTINE hpg_zps
510
511
512	SUBROUTINE hpg_sco( kt, Kmm, puu, pvv, Krhs )
513	!!---------------------------------------------------------------------
514	!! * ROUTINE hpg_sco *
515	!!
516	!! ** Method : s-coordinate case. Jacobian scheme.
517	!! The now hydrostatic pressure gradient at a given level, jk,
518	!! is computed by taking the vertical integral of the in-situ
519	!! density gradient along the model level from the suface to that
520	!! level. s-coordinates (ln_sco): a corrective term is added
521	!! to the horizontal pressure gradient :
522	!! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ]
523	!! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ]
524	!! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
525	!! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi
526	!! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj
527	!!
528	!! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
529	!!----------------------------------------------------------------------
530	INTEGER , INTENT( in ) :: kt ! ocean time-step index
531	INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
532	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
533	!!
534	INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices
535	REAL(wp) :: zcoef0, zuap, zvap, ztmp ! local scalars
536	LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables
537	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhpi, zhpj
538	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter
539	!!----------------------------------------------------------------------
540	!
541	IF( ln_wd_il ) ALLOCATE(zcpx(A2D(nn_hls)), zcpy(A2D(nn_hls)))
542	!
543	IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
544	IF( kt == nit000 ) THEN
545	IF(lwp) WRITE(numout,*)
546	IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
547	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OCE original scheme used'
548	ENDIF
549	ENDIF
550	!
551	zcoef0 = - grav * 0.5_wp
552	!
553	IF( ln_wd_il ) THEN
554	DO_2D( 0, 0, 0, 0 )
555	ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > &
556	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
557	& MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) &
558	& > rn_wdmin1 + rn_wdmin2
559	ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( &
560	& MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > &
561	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
562
563	IF(ll_tmp1) THEN
564	zcpx(ji,jj) = 1.0_wp
565	ELSE IF(ll_tmp2) THEN
566	! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here
567	zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
568	& / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) )
569	ELSE
570	zcpx(ji,jj) = 0._wp
571	END IF
572
573	ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > &
574	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
575	& MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) &
576	& > rn_wdmin1 + rn_wdmin2
577	ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( &
578	& MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > &
579	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
580
581	IF(ll_tmp1) THEN
582	zcpy(ji,jj) = 1.0_wp
583	ELSE IF(ll_tmp2) THEN
584	! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here
585	zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
586	& / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) )
587	ELSE
588	zcpy(ji,jj) = 0._wp
589	END IF
590	END_2D
591	END IF
592	!
593	DO_2D( 0, 0, 0, 0 ) ! Surface value
594	! ! hydrostatic pressure gradient along s-surfaces
595	zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) &
596	& * ( e3w(ji+1,jj ,1,Kmm) * rhd(ji+1,jj ,1) &
597	& - e3w(ji ,jj ,1,Kmm) * rhd(ji ,jj ,1) )
598	zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) &
599	& * ( e3w(ji ,jj+1,1,Kmm) * rhd(ji ,jj+1,1) &
600	& - e3w(ji ,jj ,1,Kmm) * rhd(ji ,jj ,1) )
601	! ! s-coordinate pressure gradient correction
602	zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) &
603	& * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj)
604	zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) &
605	& * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj)
606	!
607	IF( ln_wd_il ) THEN
608	zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
609	zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)
610	zuap = zuap * zcpx(ji,jj)
611	zvap = zvap * zcpy(ji,jj)
612	ENDIF
613	! ! add to the general momentum trend
614	puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1) + zuap
615	pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1) + zvap
616	END_2D
617	!
618	DO_3D( 0, 0, 0, 0, 2, jpkm1 ) ! interior value (2=<jk=<jpkm1)
619	! ! hydrostatic pressure gradient along s-surfaces
620	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) &
621	& * ( e3w(ji+1,jj,jk,Kmm) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) &
622	& - e3w(ji ,jj,jk,Kmm) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) )
623	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) &
624	& * ( e3w(ji,jj+1,jk,Kmm) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) &
625	& - e3w(ji,jj ,jk,Kmm) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) )
626	! ! s-coordinate pressure gradient correction
627	zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) ) &
628	& * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) * r1_e1u(ji,jj)
629	zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) ) &
630	& * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) * r1_e2v(ji,jj)
631	!
632	IF( ln_wd_il ) THEN
633	zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
634	zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)
635	zuap = zuap * zcpx(ji,jj)
636	zvap = zvap * zcpy(ji,jj)
637	ENDIF
638	!
639	! add to the general momentum trend
640	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk) + zuap
641	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk) + zvap
642	END_3D
643	!
644	IF( ln_wd_il ) DEALLOCATE( zcpx , zcpy )
645	!
646	END SUBROUTINE hpg_sco
647
648
649	SUBROUTINE hpg_isf( kt, Kmm, puu, pvv, Krhs )
650	!!---------------------------------------------------------------------
651	!! * ROUTINE hpg_isf *
652	!!
653	!! ** Method : s-coordinate case. Jacobian scheme.
654	!! The now hydrostatic pressure gradient at a given level, jk,
655	!! is computed by taking the vertical integral of the in-situ
656	!! density gradient along the model level from the suface to that
657	!! level. s-coordinates (ln_sco): a corrective term is added
658	!! to the horizontal pressure gradient :
659	!! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ]
660	!! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ]
661	!! add it to the general momentum trend (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)).
662	!! puu(:,:,:,Krhs) = puu(:,:,:,Krhs) - 1/e1u * zhpi
663	!! pvv(:,:,:,Krhs) = pvv(:,:,:,Krhs) - 1/e2v * zhpj
664	!! iceload is added
665	!!
666	!! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
667	!!----------------------------------------------------------------------
668	INTEGER , INTENT( in ) :: kt ! ocean time-step index
669	INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
670	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
671	!!
672	INTEGER :: ji, jj, jk ! dummy loop indices
673	INTEGER :: ikt , ikti1, iktj1 ! local integer
674	REAL(wp) :: ze3w, ze3wi1, ze3wj1 ! local scalars
675	REAL(wp) :: zcoef0, zuap, zvap ! - -
676	REAL(wp), DIMENSION(A2D(nn_hls),jpk ) :: zhpi, zhpj
677	REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: zts_top
678	REAL(wp), DIMENSION(A2D(nn_hls)) :: zrhdtop_oce
679	!!----------------------------------------------------------------------
680	!
681	zcoef0 = - grav * 0.5_wp ! Local constant initialization
682	!
683	! ! iniitialised to 0. zhpi zhpi
684	zhpi(:,:,:) = 0._wp ; zhpj(:,:,:) = 0._wp
685
686	! compute rhd at the ice/oce interface (ocean side)
687	! usefull to reduce residual current in the test case ISOMIP with no melting
688	DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
689	ikt = mikt(ji,jj)
690	zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm)
691	zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm)
692	END_2D
693	CALL eos( zts_top, risfdep, zrhdtop_oce )
694
695	! !===========================!
696	! !===== surface value =====!
697	! !===========================!
698	DO_2D( 0, 0, 0, 0 )
699	ikt = mikt(ji ,jj ) ; ze3w = e3w(ji ,jj ,ikt ,Kmm)
700	ikti1 = mikt(ji+1,jj ) ; ze3wi1 = e3w(ji+1,jj ,ikti1,Kmm)
701	iktj1 = mikt(ji ,jj+1) ; ze3wj1 = e3w(ji ,jj+1,iktj1,Kmm)
702	! ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
703	! ! we assume ISF is in isostatic equilibrium
704	zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( risfload(ji+1,jj) - risfload(ji,jj) &
705	& + 0.5_wp * ( ze3wi1 * ( rhd(ji+1,jj,ikti1) + zrhdtop_oce(ji+1,jj) ) &
706	& - ze3w * ( rhd(ji ,jj,ikt ) + zrhdtop_oce(ji ,jj) ) ) )
707	zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( risfload(ji,jj+1) - risfload(ji,jj) &
708	& + 0.5_wp * ( ze3wj1 * ( rhd(ji,jj+1,iktj1) + zrhdtop_oce(ji,jj+1) ) &
709	& - ze3w * ( rhd(ji,jj ,ikt ) + zrhdtop_oce(ji,jj ) ) ) )
710	! ! s-coordinate pressure gradient correction (=0 if z coordinate)
711	zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) &
712	& * ( gde3w(ji+1,jj,1) - gde3w(ji,jj,1) ) * r1_e1u(ji,jj)
713	zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) &
714	& * ( gde3w(ji,jj+1,1) - gde3w(ji,jj,1) ) * r1_e2v(ji,jj)
715	! ! add to the general momentum trend
716	puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1)
717	pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1)
718	END_2D
719	!
720	! !=============================!
721	! !===== interior values =====!
722	! !=============================!
723	DO_3D( 0, 0, 0, 0, 2, jpkm1 )
724	ze3w = e3w(ji ,jj ,jk,Kmm)
725	ze3wi1 = e3w(ji+1,jj ,jk,Kmm)
726	ze3wj1 = e3w(ji ,jj+1,jk,Kmm)
727	! ! hydrostatic pressure gradient along s-surfaces
728	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) &
729	& * ( ze3wi1 * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) * wmask(ji+1,jj,jk) &
730	& - ze3w * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) * wmask(ji ,jj,jk) )
731	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) &
732	& * ( ze3wj1 * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) * wmask(ji,jj+1,jk) &
733	& - ze3w * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) * wmask(ji,jj ,jk) )
734	! ! s-coordinate pressure gradient correction
735	zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) ) &
736	& * ( gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk) ) / e1u(ji,jj)
737	zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) ) &
738	& * ( gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk) ) / e2v(ji,jj)
739	! ! add to the general momentum trend
740	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
741	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk)
742	END_3D
743	!
744	END SUBROUTINE hpg_isf
745
746
747	SUBROUTINE hpg_djc( kt, Kmm, puu, pvv, Krhs )
748	!!---------------------------------------------------------------------
749	!! * ROUTINE hpg_djc *
750	!!
751	!! ** Method : Density Jacobian with Cubic polynomial scheme
752	!!
753	!! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003
754	!!----------------------------------------------------------------------
755	INTEGER , INTENT( in ) :: kt ! ocean time-step index
756	INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
757	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
758	!!
759	INTEGER :: ji, jj, jk ! dummy loop indices
760	INTEGER :: iktb, iktt ! jk indices at tracer points for top and bottom points
761	REAL(wp) :: zcoef0, zep, cffw ! temporary scalars
762	REAL(wp) :: z_grav_10, z1_12, z1_cff
763	REAL(wp) :: cffu, cffx ! " "
764	REAL(wp) :: cffv, cffy ! " "
765	LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables
766	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhpi, zhpj
767
768	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdzx, zdzy, zdzz ! Primitive grid differences ('delta_xyz')
769	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdz_i, zdz_j, zdz_k ! Harmonic average of primitive grid differences ('d_xyz')
770	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdrhox, zdrhoy, zdrhoz ! Primitive rho differences ('delta_rho')
771	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdrho_i, zdrho_j, zdrho_k ! Harmonic average of primitive rho differences ('d_rho')
772	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z_rho_i, z_rho_j, z_rho_k ! Face intergrals
773	REAL(wp), DIMENSION(A2D(nn_hls)) :: zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j ! temporary arrays
774	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter
775	!!----------------------------------------------------------------------
776	!
777	IF( ln_wd_il ) THEN
778	ALLOCATE( zcpx(A2D(nn_hls)) , zcpy(A2D(nn_hls)) )
779	DO_2D( 0, 0, 0, 0 )
780	ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > &
781	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
782	& MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) &
783	& > rn_wdmin1 + rn_wdmin2
784	ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( &
785	& MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > &
786	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
787	IF(ll_tmp1) THEN
788	zcpx(ji,jj) = 1.0_wp
789	ELSE IF(ll_tmp2) THEN
790	! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here
791	zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
792	& / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) )
793	ELSE
794	zcpx(ji,jj) = 0._wp
795	END IF
796
797	ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > &
798	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
799	& MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) &
800	& > rn_wdmin1 + rn_wdmin2
801	ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( &
802	& MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > &
803	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
804
805	IF(ll_tmp1) THEN
806	zcpy(ji,jj) = 1.0_wp
807	ELSE IF(ll_tmp2) THEN
808	! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here
809	zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
810	& / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) )
811	ELSE
812	zcpy(ji,jj) = 0._wp
813	END IF
814	END_2D
815	END IF
816
817	IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
818	IF( kt == nit000 ) THEN
819	IF(lwp) WRITE(numout,*)
820	IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend'
821	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme'
822	ENDIF
823	ENDIF
824
825	! Local constant initialization
826	zcoef0 = - grav * 0.5_wp
827	z_grav_10 = grav / 10._wp
828	z1_12 = 1.0_wp / 12._wp
829
830	!----------------------------------------------------------------------------------------
831	! 1. compute and store elementary vertical differences in provisional arrays
832	!----------------------------------------------------------------------------------------
833
834	!!bug gm Not a true bug, but... zdzz=e3w for zdzx, zdzy verify what it is really
835
836	DO_3D( 1, 1, 1, 1, 2, jpkm1 )
837	zdrhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1)
838	zdzz (ji,jj,jk) = - gde3w(ji ,jj ,jk) + gde3w(ji,jj,jk-1)
839	END_3D
840
841	!-------------------------------------------------------------------------
842	! 2. compute harmonic averages for vertical differences using eq. 5.18
843	!-------------------------------------------------------------------------
844	zep = 1.e-15
845
846	!! mb zdrho_k, zdz_k, zdrho_i, zdz_i, zdrho_j, zdz_j re-centred about the point (ji,jj,jk)
847	zdrho_k(:,:,:) = 0._wp
848	zdz_k (:,:,:) = 0._wp
849
850	DO_3D( 1, 1, 1, 1, 2, jpk-2 )
851	cffw = MAX( 2._wp * zdrhoz(ji,jj,jk) * zdrhoz(ji,jj,jk+1), 0._wp )
852	z1_cff = zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1)
853	zdrho_k(ji,jj,jk) = cffw / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
854	zdz_k(ji,jj,jk) = 2._wp * zdzz(ji,jj,jk) * zdzz(ji,jj,jk+1) &
855	& / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) )
856	END_3D
857
858	!----------------------------------------------------------------------------------
859	! 3. apply boundary conditions at top and bottom using 5.36-5.37
860	!----------------------------------------------------------------------------------
861
862	! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition
863	DO_2D( 1, 1, 1, 1 )
864	zdrho_k(ji,jj,1) = aco_bc_rhd_srf * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) - bco_bc_rhd_srf * zdrho_k(ji,jj,2)
865	zdz_k (ji,jj,1) = aco_bc_z_srf * (-gde3w(ji,jj,2) + gde3w(ji,jj,1) ) - bco_bc_z_srf * zdz_k (ji,jj,2)
866	END_2D
867
868	DO_2D( 1, 1, 1, 1 )
869	IF ( mbkt(ji,jj)>1 ) THEN
870	iktb = mbkt(ji,jj)
871	zdrho_k(ji,jj,iktb) = aco_bc_rhd_bot * ( rhd(ji,jj,iktb) - rhd(ji,jj,iktb-1) ) - bco_bc_rhd_bot * zdrho_k(ji,jj,iktb-1)
872	zdz_k (ji,jj,iktb) = aco_bc_z_bot * ( -gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_z_bot * zdz_k (ji,jj,iktb-1)
873	END IF
874	END_2D
875
876	IF ( ln_dbg_hpg ) CALL dbg_3dr( '3. zdz_k', zdz_k )
877	IF ( ln_dbg_hpg ) CALL dbg_3dr( '3. zdrho_k', zdrho_k )
878
879	!--------------------------------------------------------------
880	! 4. Compute side face integrals
881	!-------------------------------------------------------------
882
883	!! ssh replaces e3w_n ; gde3w is a depth; the formulae involve heights
884	!! rho_k stores grav * FX / rho_0
885
886	!--------------------------------------------------------------
887	! 4. a) Upper half of top-most grid box, compute and store
888	!-------------------------------------------------------------
889	! *** AY note: ssh(ji,jj,Kmm) + gde3w(ji,jj,1) = e3w(ji,jj,1,Kmm)
890	DO_2D( 0, 1, 0, 1)
891	z_rho_k(ji,jj,1) = grav * ( ssh(ji,jj,Kmm) + gde3w(ji,jj,1) ) &
892	& * ( rhd(ji,jj,1) &
893	& + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) &
894	& * ( ssh (ji,jj,Kmm) + gde3w(ji,jj,1) ) &
895	& / ( - gde3w(ji,jj,2) + gde3w(ji,jj,1) ) )
896	END_2D
897
898	!--------------------------------------------------------------
899	! 4. b) Interior faces, compute and store
900	!-------------------------------------------------------------
901
902	DO_3D( 0, 1, 0, 1, 2, jpkm1 )
903	z_rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) &
904	& * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) ) &
905	& + z_grav_10 * ( &
906	& ( zdrho_k (ji,jj,jk) - zdrho_k (ji,jj,jk-1) ) &
907	& * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) - z1_12 * ( zdz_k (ji,jj,jk) + zdz_k (ji,jj,jk-1) ) ) &
908	& - ( zdz_k (ji,jj,jk) - zdz_k (ji,jj,jk-1) ) &
909	& * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) ) &
910	& )
911	END_3D
912
913	IF ( ln_dbg_hpg ) CALL dbg_3dr( '4. z_rho_k', z_rho_k )
914
915	!----------------------------------------------------------------------------------------
916	! 5. compute and store elementary horizontal differences in provisional arrays
917	!----------------------------------------------------------------------------------------
918	zdrhox(:,:,:) = 0._wp
919	zdzx (:,:,:) = 0._wp
920	zdrhoy(:,:,:) = 0._wp
921	zdzy (:,:,:) = 0._wp
922
923	DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 )
924	zdrhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji ,jj ,jk)
925	zdzx (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji+1,jj ,jk)
926	zdrhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji ,jj ,jk)
927	zdzy (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji ,jj+1,jk)
928	END_3D
929
930	IF( nn_hls == 1 ) CALL lbc_lnk( 'dynhpg', zdrhox, 'U', -1._wp, zdzx, 'U', -1._wp, zdrhoy, 'V', -1._wp, zdzy, 'V', -1._wp )
931
932	!-------------------------------------------------------------------------
933	! 6. compute harmonic averages using eq. 5.18
934	!-------------------------------------------------------------------------
935
936	DO_3D( 0, 1, 0, 1, 1, jpkm1 )
937	cffu = MAX( 2._wp * zdrhox(ji-1,jj,jk) * zdrhox(ji,jj,jk), 0._wp )
938	z1_cff = zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk)
939	zdrho_i(ji,jj,jk) = cffu / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
940
941	cffx = MAX( 2._wp * zdzx(ji-1,jj,jk) * zdzx(ji,jj,jk), 0._wp )
942	z1_cff = zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk)
943	zdz_i(ji,jj,jk) = cffx / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
944
945	cffv = MAX( 2._wp * zdrhoy(ji,jj-1,jk) * zdrhoy(ji,jj,jk), 0._wp )
946	z1_cff = zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk)
947	zdrho_j(ji,jj,jk) = cffv / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
948
949	cffy = MAX( 2._wp * zdzy(ji,jj-1,jk) * zdzy(ji,jj,jk), 0._wp )
950	z1_cff = zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk)
951	zdz_j(ji,jj,jk) = cffy / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
952	END_3D
953
954	!!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point
955
956	!----------------------------------------------------------------------------------
957	! 6B. apply boundary conditions at side boundaries using 5.36-5.37
958	!----------------------------------------------------------------------------------
959
960	DO jk = 1, jpkm1
961	zz_drho_i(:,:) = zdrho_i(:,:,jk)
962	zz_dz_i (:,:) = zdz_i (:,:,jk)
963	zz_drho_j(:,:) = zdrho_j(:,:,jk)
964	zz_dz_j (:,:) = zdz_j (:,:,jk)
965	! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj)
966	DO_2D( 0, 0, 0, 1 )
967	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
968	zz_drho_i(ji,jj) = aco_bc_rhd_hor * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) ) - bco_bc_rhd_hor * zdrho_i(ji+1,jj,jk)
969	zz_dz_i (ji,jj) = aco_bc_z_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_z_hor * zdz_i (ji+1,jj,jk)
970	END IF
971	END_2D
972	! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj)
973	DO_2D( -1, 1, 0, 1 )
974	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
975	zz_drho_i(ji,jj) = aco_bc_rhd_hor * ( rhd (ji,jj,jk) - rhd (ji-1,jj,jk) ) - bco_bc_rhd_hor * zdrho_i(ji-1,jj,jk)
976	zz_dz_i (ji,jj) = aco_bc_z_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_z_hor * zdz_i (ji-1,jj,jk)
977	END IF
978	END_2D
979	! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi)
980	DO_2D( 0, 1, 0, 0 )
981	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
982	zz_drho_j(ji,jj) = aco_bc_rhd_hor * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) ) - bco_bc_rhd_hor * zdrho_j(ji,jj+1,jk)
983	zz_dz_j (ji,jj) = aco_bc_z_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_z_hor * zdz_j (ji,jj+1,jk)
984	END IF
985	END_2D
986	! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi)
987	DO_2D( 0, 1, -1, 1 )
988	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
989	zz_drho_j(ji,jj) = aco_bc_rhd_hor * ( rhd (ji,jj,jk) - rhd (ji,jj-1,jk) ) - bco_bc_rhd_hor * zdrho_j(ji,jj-1,jk)
990	zz_dz_j (ji,jj) = aco_bc_z_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_z_hor * zdz_j (ji,jj-1,jk)
991	END IF
992	END_2D
993	zdrho_i(:,:,jk) = zz_drho_i(:,:)
994	zdz_i (:,:,jk) = zz_dz_i (:,:)
995	zdrho_j(:,:,jk) = zz_drho_j(:,:)
996	zdz_j (:,:,jk) = zz_dz_j (:,:)
997	END DO ! jk
998
999	IF ( ln_dbg_hpg ) THEN
1000	CALL dbg_3dr( '6. zdrho_i', zdrho_i )
1001	CALL dbg_3dr( '6. zdrho_j', zdrho_j )
1002	END IF
1003
1004	!--------------------------------------------------------------
1005	! 7. Calculate integrals on upper/lower faces
1006	!-------------------------------------------------------------
1007
1008	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
1009	! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height)
1010	z_rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) &
1011	& * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )
1012	IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN
1013	z_rho_i(ji,jj,jk) = z_rho_i(ji,jj,jk) - z_grav_10 * ( &
1014	& ( zdrho_i (ji+1,jj,jk) - zdrho_i (ji,jj,jk) ) &
1015	& * ( - gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_i (ji+1,jj,jk) + zdz_i (ji,jj,jk) ) ) &
1016	& - ( zdz_i (ji+1,jj,jk) - zdz_i (ji,jj,jk) ) &
1017	& * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) ) &
1018	& )
1019	END IF
1020
1021	z_rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) &
1022	& * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )
1023	IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN
1024	z_rho_j(ji,jj,jk) = z_rho_j(ji,jj,jk) - z_grav_10 * ( &
1025	& ( zdrho_j (ji,jj+1,jk) - zdrho_j (ji,jj,jk) ) &
1026	& * ( - gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_j (ji,jj+1,jk) + zdz_j (ji,jj,jk) ) ) &
1027	& - ( zdz_j (ji,jj+1,jk) - zdz_j (ji,jj,jk) ) &
1028	& * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) ) &
1029	& )
1030	END IF
1031	END_3D
1032
1033	IF ( ln_dbg_hpg ) THEN
1034	CALL dbg_3dr( '7. z_rho_i', z_rho_i )
1035	CALL dbg_3dr( '7. z_rho_j', z_rho_j )
1036	END IF
1037
1038	!--------------------------------------------------------------
1039	! 8. Integrate in the vertical
1040	!-------------------------------------------------------------
1041	!
1042	! ---------------
1043	! Surface value
1044	! ---------------
1045	DO_2D( 0, 0, 0, 0 )
1046	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)
1047	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)
1048	IF( ln_wd_il ) THEN
1049	zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
1050	zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)
1051	ENDIF
1052	! add to the general momentum trend
1053	puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) + zhpi(ji,jj,1)
1054	pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) + zhpj(ji,jj,1)
1055	END_2D
1056
1057	! ----------------
1058	! interior value (2=<jk=<jpkm1)
1059	! ----------------
1060	DO_3D( 0, 0, 0, 0, 2, jpkm1 )
1061	! hydrostatic pressure gradient along s-surfaces
1062	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &
1063	& + ( ( z_rho_k(ji,jj,jk) - z_rho_k(ji+1,jj,jk ) ) &
1064	& - ( z_rho_i(ji,jj,jk) - z_rho_i(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj)
1065	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &
1066	& + ( ( z_rho_k(ji,jj,jk) - z_rho_k(ji,jj+1,jk ) ) &
1067	& -( z_rho_j(ji,jj,jk) - z_rho_j(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj)
1068	IF( ln_wd_il ) THEN
1069	zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
1070	zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)
1071	ENDIF
1072	! add to the general momentum trend
1073	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj,jk)
1074	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj,jk)
1075	END_3D
1076	!
1077
1078	IF ( ln_dbg_hpg ) THEN
1079	CALL dbg_3dr( '8. zhpi', zhpi )
1080	CALL dbg_3dr( '8. zhpj', zhpj )
1081	END IF
1082
1083	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
1084	!
1085	END SUBROUTINE hpg_djc
1086
1087
1088	SUBROUTINE hpg_prj( kt, Kmm, puu, pvv, Krhs )
1089	!!---------------------------------------------------------------------
1090	!! * ROUTINE hpg_prj *
1091	!!
1092	!! ** Method : s-coordinate case.
1093	!! A Pressure-Jacobian horizontal pressure gradient method
1094	!! based on the constrained cubic-spline interpolation for
1095	!! all vertical coordinate systems
1096	!!
1097	!! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now hydrastatic pressure trend
1098	!!----------------------------------------------------------------------
1099	INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear
1100	INTEGER , INTENT( in ) :: kt ! ocean time-step index
1101	INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
1102	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
1103	!!
1104	INTEGER :: ji, jj, jk, jkk ! dummy loop indices
1105	REAL(wp) :: zcoef0, znad ! local scalars
1106	!
1107	!! The local variables for the correction term
1108	INTEGER :: jk1, jis, jid, jjs, jjd
1109	LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables
1110	REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps
1111	REAL(wp) :: zrhdt1
1112	REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2
1113	REAL(wp), DIMENSION(A2D(nn_hls)) :: zpgu, zpgv ! 2D workspace
1114	REAL(wp), DIMENSION(A2D(nn_hls)) :: zsshu_n, zsshv_n
1115	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdept, zrhh
1116	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
1117	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter
1118	!!----------------------------------------------------------------------
1119	!
1120	IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
1121	IF( kt == nit000 ) THEN
1122	IF(lwp) WRITE(numout,*)
1123	IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend'
1124	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, cubic spline pressure Jacobian'
1125	ENDIF
1126	ENDIF
1127
1128	! Local constant initialization
1129	zcoef0 = - grav
1130	znad = 1._wp
1131	IF( ln_linssh ) znad = 1._wp
1132	!
1133	! ---------------
1134	! Surface pressure gradient to be removed
1135	! ---------------
1136	DO_2D( 0, 0, 0, 0 )
1137	zpgu(ji,jj) = - grav * ( ssh(ji+1,jj,Kmm) - ssh(ji,jj,Kmm) ) * r1_e1u(ji,jj)
1138	zpgv(ji,jj) = - grav * ( ssh(ji,jj+1,Kmm) - ssh(ji,jj,Kmm) ) * r1_e2v(ji,jj)
1139	END_2D
1140	!
1141	IF( ln_wd_il ) THEN
1142	ALLOCATE( zcpx(A2D(nn_hls)) , zcpy(A2D(nn_hls)) )
1143	DO_2D( 0, 0, 0, 0 )
1144	ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > &
1145	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
1146	& MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) > &
1147	& rn_wdmin1 + rn_wdmin2
1148	ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. &
1149	& ( MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > &
1150	& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 )
1151
1152	IF(ll_tmp1) THEN
1153	zcpx(ji,jj) = 1.0_wp
1154	ELSE IF(ll_tmp2) THEN
1155	! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here
1156	zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
1157	& / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) )
1158	zcpx(ji,jj) = MAX(MIN( zcpx(ji,jj) , 1.0_wp),0.0_wp)
1159	ELSE
1160	zcpx(ji,jj) = 0._wp
1161	END IF
1162
1163	ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > &
1164	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. &
1165	& MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) > &
1166	& rn_wdmin1 + rn_wdmin2
1167	ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. &
1168	& ( MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > &
1169	& MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 )
1170
1171	IF(ll_tmp1) THEN
1172	zcpy(ji,jj) = 1.0_wp
1173	ELSE IF(ll_tmp2) THEN
1174	! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here
1175	zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) &
1176	& / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) )
1177	zcpy(ji,jj) = MAX(MIN( zcpy(ji,jj) , 1.0_wp),0.0_wp)
1178	ELSE
1179	zcpy(ji,jj) = 0._wp
1180	ENDIF
1181	END_2D
1182	ENDIF
1183
1184	! Clean 3-D work arrays
1185	zhpi(:,:,:) = 0._wp
1186	zrhh(:,:,:) = rhd(A2D(nn_hls),:)
1187
1188	! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate
1189	DO_2D( 1, 1, 1, 1 )
1190	jk = mbkt(ji,jj)
1191	IF( jk <= 1 ) THEN ; zrhh(ji,jj, : ) = 0._wp
1192	ELSEIF( jk == 2 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk)
1193	ELSEIF( jk < jpkm1 ) THEN
1194	DO jkk = jk+1, jpk
1195	zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk ), gde3w(ji,jj,jkk-1), &
1196	& gde3w(ji,jj,jkk-2), zrhh (ji,jj,jkk-1), zrhh(ji,jj,jkk-2))
1197	END DO
1198	ENDIF
1199	END_2D
1200
1201	! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)"
1202	DO_2D( 1, 1, 1, 1 )
1203	zdept(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) - ssh(ji,jj,Kmm)
1204	END_2D
1205
1206	DO_3D( 1, 1, 1, 1, 2, jpk )
1207	zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + e3w(ji,jj,jk,Kmm)
1208	END_3D
1209
1210	fsp(:,:,:) = zrhh (:,:,:)
1211	xsp(:,:,:) = zdept(:,:,:)
1212
1213	! Construct the vertical density profile with the
1214	! constrained cubic spline interpolation
1215	! rho(z) = asp + bspz + cspz^2 + dsp*z^3
1216	CALL cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type )
1217
1218	! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)"
1219	DO_2D( 0, 1, 0, 1 )
1220	zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), &
1221	& csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm)
1222
1223	! assuming linear profile across the top half surface layer
1224	zhpi(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1
1225	END_2D
1226
1227	! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)"
1228	DO_3D( 0, 1, 0, 1, 2, jpkm1 )
1229	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + &
1230	& integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), &
1231	& asp (ji,jj,jk-1), bsp (ji,jj,jk-1), &
1232	& csp (ji,jj,jk-1), dsp (ji,jj,jk-1) )
1233	END_3D
1234
1235	! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1)
1236
1237	! Prepare zsshu_n and zsshv_n
1238	DO_2D( 0, 0, 0, 0 )
1239	!!gm BUG ? if it is ssh at u- & v-point then it should be:
1240	! zsshu_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji+1,jj) * ssh(ji+1,jj,Kmm)) * &
1241	! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp
1242	! zsshv_n(ji,jj) = (e1e2t(ji,jj) * ssh(ji,jj,Kmm) + e1e2t(ji,jj+1) * ssh(ji,jj+1,Kmm)) * &
1243	! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp
1244	!!gm not this:
1245	zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * &
1246	& r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp
1247	zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * &
1248	& r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp
1249	END_2D
1250
1251	DO_2D( 0, 0, 0, 0 )
1252	zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) )
1253	zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) )
1254	END_2D
1255
1256	DO_3D( 0, 0, 0, 0, 2, jpkm1 )
1257	zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm)
1258	zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm)
1259	END_3D
1260
1261	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
1262	zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm)
1263	zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm)
1264	END_3D
1265
1266	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
1267	zu(ji,jj,jk) = MIN( zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) )
1268	zu(ji,jj,jk) = MAX( zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) )
1269	zv(ji,jj,jk) = MIN( zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) )
1270	zv(ji,jj,jk) = MAX( zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) )
1271	END_3D
1272
1273
1274	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
1275	zpwes = 0._wp; zpwed = 0._wp
1276	zpnss = 0._wp; zpnsd = 0._wp
1277	zuijk = zu(ji,jj,jk)
1278	zvijk = zv(ji,jj,jk)
1279
1280	!!!!! for u equation
1281	IF( jk <= mbku(ji,jj) ) THEN
1282	IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN
1283	jis = ji + 1; jid = ji
1284	ELSE
1285	jis = ji; jid = ji +1
1286	ENDIF
1287
1288	! integrate the pressure on the shallow side
1289	jk1 = jk
1290	DO WHILE ( -zdept(jis,jj,jk1) > zuijk )
1291	IF( jk1 == mbku(ji,jj) ) THEN
1292	zuijk = -zdept(jis,jj,jk1)
1293	EXIT
1294	ENDIF
1295	zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk)
1296	zpwes = zpwes + &
1297	integ_spline(zdept(jis,jj,jk1), zdeps, &
1298	asp(jis,jj,jk1), bsp(jis,jj,jk1), &
1299	csp(jis,jj,jk1), dsp(jis,jj,jk1))
1300	jk1 = jk1 + 1
1301	END DO
1302
1303	! integrate the pressure on the deep side
1304	jk1 = jk
1305	DO WHILE ( -zdept(jid,jj,jk1) < zuijk )
1306	IF( jk1 == 1 ) THEN
1307	zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad)
1308	zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), &
1309	bsp(jid,jj,1) , csp(jid,jj,1), &
1310	dsp(jid,jj,1)) * zdeps
1311	zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps
1312	EXIT
1313	ENDIF
1314	zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk)
1315	zpwed = zpwed + &
1316	integ_spline(zdeps, zdept(jid,jj,jk1), &
1317	asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), &
1318	csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) )
1319	jk1 = jk1 - 1
1320	END DO
1321
1322	! update the momentum trends in u direction
1323	zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) )
1324	IF( .NOT.ln_linssh ) THEN
1325	zdpdx2 = zcoef0 * r1_e1u(ji,jj) * &
1326	& ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) )
1327	ELSE
1328	zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)
1329	ENDIF
1330	IF( ln_wd_il ) THEN
1331	zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj)
1332	zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj)
1333	ENDIF
1334	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk)
1335	ENDIF
1336
1337	!!!!! for v equation
1338	IF( jk <= mbkv(ji,jj) ) THEN
1339	IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN
1340	jjs = jj + 1; jjd = jj
1341	ELSE
1342	jjs = jj ; jjd = jj + 1
1343	ENDIF
1344
1345	! integrate the pressure on the shallow side
1346	jk1 = jk
1347	DO WHILE ( -zdept(ji,jjs,jk1) > zvijk )
1348	IF( jk1 == mbkv(ji,jj) ) THEN
1349	zvijk = -zdept(ji,jjs,jk1)
1350	EXIT
1351	ENDIF
1352	zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk)
1353	zpnss = zpnss + &
1354	integ_spline(zdept(ji,jjs,jk1), zdeps, &
1355	asp(ji,jjs,jk1), bsp(ji,jjs,jk1), &
1356	csp(ji,jjs,jk1), dsp(ji,jjs,jk1) )
1357	jk1 = jk1 + 1
1358	END DO
1359
1360	! integrate the pressure on the deep side
1361	jk1 = jk
1362	DO WHILE ( -zdept(ji,jjd,jk1) < zvijk )
1363	IF( jk1 == 1 ) THEN
1364	zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad)
1365	zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), &
1366	bsp(ji,jjd,1) , csp(ji,jjd,1), &
1367	dsp(ji,jjd,1) ) * zdeps
1368	zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps
1369	EXIT
1370	ENDIF
1371	zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk)
1372	zpnsd = zpnsd + &
1373	integ_spline(zdeps, zdept(ji,jjd,jk1), &
1374	asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), &
1375	csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
1376	jk1 = jk1 - 1
1377	END DO
1378
1379	! update the momentum trends in v direction
1380	zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) )
1381	IF( .NOT.ln_linssh ) THEN
1382	zdpdy2 = zcoef0 * r1_e2v(ji,jj) * &
1383	( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) )
1384	ELSE
1385	zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )
1386	ENDIF
1387	IF( ln_wd_il ) THEN
1388	zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj)
1389	zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj)
1390	ENDIF
1391
1392	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk)
1393	ENDIF
1394	!
1395	END_3D
1396	!
1397	IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
1398	!
1399	END SUBROUTINE hpg_prj
1400
1401
1402	SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type )
1403	!!----------------------------------------------------------------------
1404	!! * ROUTINE cspline *
1405	!!
1406	!! ** Purpose : constrained cubic spline interpolation
1407	!!
1408	!! ** Method : f(x) = asp + bspx + cspx^2 + dsp*x^3
1409	!!
1410	!! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation
1411	!!----------------------------------------------------------------------
1412	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: fsp, xsp ! value and coordinate
1413	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT( out) :: asp, bsp, csp, dsp ! coefficients of the interpoated function
1414	INTEGER , INTENT(in ) :: polynomial_type ! 1: cubic spline ; 2: Linear
1415	!
1416	INTEGER :: ji, jj, jk ! dummy loop indices
1417	REAL(wp) :: zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp
1418	REAL(wp) :: zdxtmp1, zdxtmp2, zalpha
1419	REAL(wp) :: zdf(jpk)
1420	!!----------------------------------------------------------------------
1421	!
1422	IF (polynomial_type == 1) THEN ! Constrained Cubic Spline
1423	DO_2D( 1, 1, 1, 1 )
1424	!!Fritsch&Butland's method, 1984 (preferred, but more computation)
1425	! DO jk = 2, jpkm1-1
1426	! zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1)
1427	! zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1428	! zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1
1429	! zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2
1430	!
1431	! zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp
1432	!
1433	! IF(zdf1 * zdf2 <= 0._wp) THEN
1434	! zdf(jk) = 0._wp
1435	! ELSE
1436	! zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 )
1437	! ENDIF
1438	! END DO
1439
1440	!!Simply geometric average
1441	DO jk = 2, jpk-2
1442	zdf1 = (fsp(ji,jj,jk ) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk ) - xsp(ji,jj,jk-1))
1443	zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk )) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk ))
1444
1445	IF(zdf1 * zdf2 <= 0._wp) THEN
1446	zdf(jk) = 0._wp
1447	ELSE
1448	zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2)
1449	ENDIF
1450	END DO
1451
1452	zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / &
1453	& ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2)
1454	zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / &
1455	& ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpk - 2)
1456
1457	DO jk = 1, jpk-2
1458	zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1459	ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp
1460	ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp
1461	zddf1 = -2._wp * ztmp1 + ztmp2
1462	ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp
1463	zddf2 = 2._wp * ztmp1 - ztmp2
1464
1465	dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp
1466	csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp
1467	bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - &
1468	& csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - &
1469	& dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - &
1470	& xsp(ji,jj,jk+1) * xsp(ji,jj,jk))
1471	asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + &
1472	& (xsp(ji,jj,jk) * (csp(ji,jj,jk) + &
1473	& dsp(ji,jj,jk) * xsp(ji,jj,jk))))
1474	END DO
1475	END_2D
1476
1477	ELSEIF ( polynomial_type == 2 ) THEN ! Linear
1478	DO_3D( 1, 1, 1, 1, 1, jpk-2 )
1479	zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1480	ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk)
1481
1482	dsp(ji,jj,jk) = 0._wp
1483	csp(ji,jj,jk) = 0._wp
1484	bsp(ji,jj,jk) = ztmp1 / zdxtmp
1485	asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk)
1486	END_3D
1487	!
1488	ELSE
1489	CALL ctl_stop( 'invalid polynomial type in cspline' )
1490	ENDIF
1491	!
1492	END SUBROUTINE cspline
1493
1494
1495	FUNCTION interp1(x, xl, xr, fl, fr) RESULT(f)
1496	!!----------------------------------------------------------------------
1497	!! * ROUTINE interp1 *
1498	!!
1499	!! ** Purpose : 1-d linear interpolation
1500	!!
1501	!! ** Method : interpolation is straight forward
1502	!! extrapolation is also permitted (no value limit)
1503	!!----------------------------------------------------------------------
1504	REAL(wp), INTENT(in) :: x, xl, xr, fl, fr
1505	REAL(wp) :: f ! result of the interpolation (extrapolation)
1506	REAL(wp) :: zdeltx
1507	!!----------------------------------------------------------------------
1508	!
1509	zdeltx = xr - xl
1510	IF( abs(zdeltx) <= 10._wp * EPSILON(x) ) THEN
1511	f = 0.5_wp * (fl + fr)
1512	ELSE
1513	f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx
1514	ENDIF
1515	!
1516	END FUNCTION interp1
1517
1518
1519	FUNCTION interp2( x, a, b, c, d ) RESULT(f)
1520	!!----------------------------------------------------------------------
1521	!! * ROUTINE interp1 *
1522	!!
1523	!! ** Purpose : 1-d constrained cubic spline interpolation
1524	!!
1525	!! ** Method : cubic spline interpolation
1526	!!
1527	!!----------------------------------------------------------------------
1528	REAL(wp), INTENT(in) :: x, a, b, c, d
1529	REAL(wp) :: f ! value from the interpolation
1530	!!----------------------------------------------------------------------
1531	!
1532	f = a + x* ( b + x * ( c + d * x ) )
1533	!
1534	END FUNCTION interp2
1535
1536
1537	FUNCTION interp3( x, a, b, c, d ) RESULT(f)
1538	!!----------------------------------------------------------------------
1539	!! * ROUTINE interp1 *
1540	!!
1541	!! ** Purpose : Calculate the first order of derivative of
1542	!! a cubic spline function y=a+bx+cx^2+d*x^3
1543	!!
1544	!! ** Method : f=dy/dx=b+2cx+3dx^2
1545	!!
1546	!!----------------------------------------------------------------------
1547	REAL(wp), INTENT(in) :: x, a, b, c, d
1548	REAL(wp) :: f ! value from the interpolation
1549	!!----------------------------------------------------------------------
1550	!
1551	f = b + x * ( 2._wp * c + 3._wp * d * x)
1552	!
1553	END FUNCTION interp3
1554
1555
1556	FUNCTION integ_spline( xl, xr, a, b, c, d ) RESULT(f)
1557	!!----------------------------------------------------------------------
1558	!! * ROUTINE interp1 *
1559	!!
1560	!! ** Purpose : 1-d constrained cubic spline integration
1561	!!
1562	!! ** Method : integrate polynomial a+bx+cx^2+dx^3 from xl to xr
1563	!!
1564	!!----------------------------------------------------------------------
1565	REAL(wp), INTENT(in) :: xl, xr, a, b, c, d
1566	REAL(wp) :: za1, za2, za3
1567	REAL(wp) :: f ! integration result
1568	!!----------------------------------------------------------------------
1569	!
1570	za1 = 0.5_wp * b
1571	za2 = c / 3.0_wp
1572	za3 = 0.25_wp * d
1573	!
1574	f = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - &
1575	& xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) )
1576	!
1577	END FUNCTION integ_spline
1578
1579	SUBROUTINE hpg_djr( kt, Kmm, puu, pvv, Krhs )
1580	!!---------------------------------------------------------------------
1581	!! * ROUTINE hpg_djr *
1582	!!
1583	!! ** Method : Density Jacobian with Cubic polynomial scheme subtracting a local reference profile (pmr is profile minus reference)
1584	!! This code assumes a 2-point halo
1585	!!
1586	!! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003
1587	!!----------------------------------------------------------------------
1588	INTEGER , INTENT( in ) :: kt ! ocean time-step index
1589	INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
1590	REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
1591	!!
1592	INTEGER :: ji, jj, jk, jr ! loop indices
1593	INTEGER :: jn_hor_pts ! number of points in the horizontal stencil
1594	INTEGER :: j_uv ! 1 for u-cell; 2 for v-cell
1595	INTEGER :: iktb, iktt ! jk indices at tracer points for top and bottom points
1596	INTEGER :: jia, jib, jja, jjb !
1597	INTEGER :: jir, jjr ! reference (expand)
1598	INTEGER :: jio, jjo ! offset (expand)
1599
1600	REAL(wp) :: z_grav_10, z1_12 ! constants
1601	REAL(wp) :: zhta, zhtb ! temporary scalars
1602	REAL(wp) :: zcoef0, zep, cffw ! " "
1603	REAL(wp) :: aco, bco ! " "
1604	REAL(wp) :: cffu, cffx, z1_cff ! " "
1605	REAL(wp) :: cffv, cffy ! " "
1606	REAL(wp) :: cff_31, cff_42 ! " "
1607	LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables
1608
1609	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: ztmp, zdz_i, zdz_j, zdz_k ! Harmonic average of primitive grid differences
1610	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zrhd_ref ! Reference density
1611	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zz_ref ! Reference heights
1612	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdrhd_k_ref ! Harmonic average of primitive differences for reference field
1613	REAL(wp), DIMENSION(A2D(nn_hls),jpk,4) :: z_rhd_pmr ! rhd_prm = rhd - rhd_ref (values on the original grid)
1614	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdrhd_k_pmr !
1615	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z_lmr_k ! left minus right density integrals on vertical faces
1616
1617	INTEGER, DIMENSION(A2D(nn_hls)) :: jk_bot_ref ! bottom levels in the reference profile
1618	REAL(wp), DIMENSION(A2D(nn_hls)) :: zdzx, zdzy ! primitive differences in x and y
1619	REAL(wp), DIMENSION(A2D(nn_hls)) :: zz_dz_i, zz_dz_j
1620	REAL(wp), DIMENSION(A2D(nn_hls)) :: zdrhd_21, zdrhd_32, zdrhd_43
1621	REAL(wp), DIMENSION(A2D(nn_hls),2) :: zz_drhd_ij, zdrhd_ij
1622	REAL(wp), DIMENSION(A2D(nn_hls)) :: z_low_ij, z_upp_ij
1623	REAL(wp), DIMENSION(A2D(nn_hls)) :: zhpi, zhpj
1624
1625	!!----------------------------------------------------------------------
1626
1627	IF( kt == nit000 ) THEN
1628	IF(lwp) WRITE(numout,*)
1629	IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend'
1630	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme'
1631	ENDIF
1632
1633	! Local constant initialization
1634	zcoef0 = - grav * 0.5_wp
1635	z_grav_10 = grav / 10._wp
1636	z1_12 = 1.0_wp / 12._wp
1637	zep = 1.e-15
1638
1639	jn_hor_pts = 4 ! 4 points in the horizontal stencil
1640
1641	!------------------------------------------------------------------------------------------------------
1642	! 1. calculate harmonic averages of differences for z (grid heights) in i, j and k directions
1643	!------------------------------------------------------------------------------------------------------
1644
1645	!------------------------------------------------------------------------------------------------------
1646	! 1.1 compute and store elementary vertical differences then harmonic averages for z using eqn 5.18
1647	! Full domain covered so that _ref profiles can be taken from zdz_k
1648	!------------------------------------------------------------------------------------------------------
1649
1650	DO_3D( 2, 2, 2, 2, 2, jpk )
1651	ztmp (ji,jj,jk) = - gde3w(ji ,jj ,jk) + gde3w(ji,jj,jk-1)
1652	END_3D
1653
1654	zdz_k (:,:,1) = 0._wp ! jk index changed from : to 1 to make computationally less wasteful
1655
1656	DO_3D( 2, 2, 2, 2, 2, jpk-1 )
1657	zdz_k(ji,jj,jk) = 2._wp * ztmp(ji,jj,jk) * ztmp(ji,jj,jk+1) &
1658	& / ( ztmp(ji,jj,jk) + ztmp(ji,jj,jk+1) )
1659	END_3D
1660
1661	!----------------------------------------------------------------------------------
1662	! 1.2 apply boundary conditions at top and bottom using 5.36-5.37
1663	!----------------------------------------------------------------------------------
1664
1665	! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition
1666	zdz_k (:,:,1) = aco_bc_z_srf * (-gde3w(:,:,2) + gde3w(:,:,1) ) - bco_bc_z_srf * zdz_k (:,:,2)
1667
1668	DO_2D( 2, 2, 2, 2 )
1669	iktb = mbkt(ji,jj)
1670	IF ( iktb > 1 ) THEN
1671	zdz_k (ji,jj,iktb) = aco_bc_z_bot * (-gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_z_bot * zdz_k (ji,jj,iktb-1)
1672	END IF
1673	END_2D
1674
1675	!----------------------------------------------------------------------------------------
1676	! 1.3 compute and store elementary horizontal differences then harmonic averages for z using eqn 5.18
1677	!----------------------------------------------------------------------------------------
1678
1679	DO jk = 1, jpkm1
1680	DO_2D( 1, 1, 1, 1 )
1681	zdzx (ji,jj) = - gde3w(ji+1,jj ,jk) + gde3w(ji,jj,jk )
1682	zdzy (ji,jj) = - gde3w(ji ,jj+1,jk) + gde3w(ji,jj,jk )
1683	END_2D
1684
1685	! this subroutine requires a 2-point halo CALL lbc_lnk_multi( 'dynhpg', zdzx, 'U', 1., zdzy, 'V', 1. )
1686
1687	DO_2D( 0, 1, 0, 1 )
1688	cffx = MAX( 2._wp * zdzx (ji-1,jj) * zdzx (ji,jj), 0._wp)
1689	z1_cff = zdzx(ji-1,jj) + zdzx(ji,jj)
1690	zdz_i(ji,jj,jk) = cffx / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
1691
1692	cffy = 2._wp * zdzy (ji ,jj-1) * zdzy (ji,jj )
1693	z1_cff = zdzy(ji,jj-1) + zdzy(ji,jj)
1694	zdz_j(ji,jj,jk) = cffy / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
1695	END_2D
1696
1697	!----------------------------------------------------------------------------------
1698	! 1.4 apply boundary conditions at sides using 5.36-5.37
1699	!----------------------------------------------------------------------------------
1700
1701	DO_2D( 0, 1, 0, 1)
1702	! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj)
1703	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
1704	zdz_i(ji,jj,jk) = aco_bc_z_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_z_hor * zz_dz_i(ji+1,jj)
1705	END IF
1706	! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj)
1707	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
1708	zdz_i(ji,jj,jk) = aco_bc_z_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_z_hor * zz_dz_i(ji-1,jj)
1709	END IF
1710	! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi)
1711	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
1712	zdz_j(ji,jj,jk) = aco_bc_z_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_z_hor * zz_dz_j(ji,jj+1)
1713	END IF
1714	! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi)
1715	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
1716	zdz_j(ji,jj,jk) = aco_bc_z_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_z_hor * zz_dz_j(ji,jj-1)
1717	END IF
1718	END_2D
1719	END DO ! k
1720
1721	IF ( ln_dbg_hpg ) THEN
1722	CALL dbg_3dr( '1.4 gde3w', gde3w )
1723	CALL dbg_3dr( '1.4 zdz_i', zdz_i )
1724	CALL dbg_3dr( '1.4 zdz_j', zdz_j )
1725	CALL dbg_3dr( '1.4 zdz_k', zdz_k )
1726	CALL dbg_2di( 'mbkt', mbkt)
1727	END IF
1728	!----------------------------------------------------------------------------------------
1729	! 2. Start loop over the u and v components and find the reference profile
1730	! The loop ends in section 5.4
1731	!----------------------------------------------------------------------------------------
1732
1733	DO j_uv = 1, 2 ! j_uv = 1 is for u-cell ; j_uv = 2 for v-cell
1734
1735	!----------------------------------------------------------------------------------------
1736	! 2.1 find reference profiles zrhd_ref and zz_ref and the bottom level of the reference profile
1737	!----------------------------------------------------------------------------------------
1738
1739	IF ( ln_dbg_hpg ) CALL dbg_2dr( '2.1 ht_0', ht_0 )
1740	CALL calc_rhd_ref(j_uv, jn_hor_pts, zrhd_ref, zz_ref, jk_bot_ref) ! Uses rhd (IN) to calculate all other fields (OUT)
1741
1742	IF ( ln_dbg_hpg ) THEN
1743	CALL dbg_3dr( '2.1 rhd', rhd )
1744	CALL dbg_3dr( '2.1 zrhd_ref', zrhd_ref )
1745	CALL dbg_3dr( '2.1 zz_ref', zz_ref )
1746	CALL dbg_2di( '2.1 jk_bot_ref', jk_bot_ref )
1747	END IF
1748
1749	!--------------------------------------------------------------------------------------------------------
1750	! 2.2 IF ln_hpg_djr_ref_ccs compute zdrhd_k_ref then set bcs at top & bottom
1751	! (bcs not needed for simple cubic off-centred at boundaries)
1752	!--------------------------------------------------------------------------------------------------------
1753
1754	IF ( ln_hpg_djr_ref_ccs ) THEN
1755	CALL calc_drhd_k(zrhd_ref, jk_bot_ref, zdrhd_k_ref)
1756	IF ( ln_dbg_hpg ) CALL dbg_3dr( '2.3 zdrhd_k_ref', zdrhd_k_ref )
1757	END IF ! ln_hpg_djr_ref_ccs
1758
1759	!----------------------------------------------------------------------------------------
1760	! 3. interpolate reference profiles to target profiles and form difference profiles z_rhd_pmr
1761	!----------------------------------------------------------------------------------------
1762
1763	DO jr = 1, 4
1764	IF ( j_uv == 1 ) THEN
1765	jio = jr - 2 ! range of jio is -1 to 2
1766	jjo = 0
1767	ELSE
1768	jio = 0
1769	jjo = jr - 2
1770	END IF
1771
1772	IF ( ln_hpg_djr_ref_ccs ) THEN
1773	CALL ref_to_tgt_ccs ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )
1774	ELSE
1775	CALL ref_to_tgt_cub ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )
1776	END IF
1777
1778	IF ( ln_dbg_hpg ) CALL dbg_3dr( '3. z_rhd_pmr', z_rhd_pmr(:,:,:,jr) )
1779
1780	END DO
1781
1782	!----------------------------------------------------------------------------------------
1783	! 4. Calculations for side-face integrals
1784	!----------------------------------------------------------------------------------------
1785
1786	!----------------------------------------------------------------------------------------
1787	! 4.1 compute and store elementary vertical differences then harmonic averages
1788	! based on z_rhd_pmr arrays (zdz_k has already been calculated)
1789	!----------------------------------------------------------------------------------------
1790
1791	! start loop over the two side faces jr = 2 "left" face; jr = 3 "right" face
1792	DO jr = 2, 3
1793
1794	IF ( j_uv == 1 ) THEN
1795	jio = jr - 2
1796	jjo = 0
1797	ELSE
1798	jio = 0
1799	jjo = jr - 2
1800	END IF
1801
1802	DO_3D( 0, 0, 0, 0, 2, jpk )
1803	ztmp(ji,jj,jk) = z_rhd_pmr (ji,jj,jk,jr) - z_rhd_pmr(ji,jj,jk-1,jr)
1804	END_3D
1805
1806	zdrhd_k_pmr(:,:,:) = 0._wp ! should be unnecessary
1807
1808	DO_3D( 0, 0, 0, 0, 2, jpk-1 )
1809	cffw = MAX( 2._wp * ztmp(ji,jj,jk) * ztmp(ji,jj,jk+1), 0._wp )
1810	z1_cff = ztmp(ji,jj,jk) + ztmp(ji,jj,jk+1)
1811	zdrhd_k_pmr(ji,jj,jk) = cffw / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
1812	END_3D
1813
1814	! apply boundary conditions at top and bottom
1815	DO_2D( 0, 0, 0, 0 )
1816	zdrhd_k_pmr(ji,jj,1) = aco_bc_rhd_srf * ( z_rhd_pmr(ji,jj,2,jr) - z_rhd_pmr(ji,jj,1,jr) ) - bco_bc_rhd_srf * zdrhd_k_pmr(ji,jj,2)
1817	iktb = mbkt(ji+jio,jj+jjo)
1818	IF ( iktb > 1 ) THEN
1819	zdrhd_k_pmr(ji,jj,iktb) = aco_bc_rhd_bot * (z_rhd_pmr(ji,jj,iktb,jr) - z_rhd_pmr(ji,jj,iktb-1,jr) ) - bco_bc_rhd_bot * zdrhd_k_pmr(ji,jj,iktb-1)
1820	END IF
1821	END_2D
1822
1823
1824	!--------------------------------------------------------------
1825	! 4.2 Upper half of top-most grid box, compute and store
1826	!-------------------------------------------------------------
1827	!! ssh replaces e3w_n ; gde3w is a depth; the formulae involve heights
1828	!! rho_k stores grav * FX / rho_0
1829	!! *** AY note: ssh(ji,jj,Kmm) + gde3w(ji,jj,1) = e3w(ji,jj,1,Kmm)
1830	DO_2D( 0, 0, 0, 0)
1831	ztmp(ji,jj,1) = grav * ( ssh(ji+jio,jj+jjo,Kmm) + gde3w(ji+jio,jj+jjo,1) ) &
1832	& * ( z_rhd_pmr(ji,jj,1,jr) &
1833	& + 0.5_wp * ( z_rhd_pmr(ji,jj,2,jr) - z_rhd_pmr(ji,jj,1,jr) ) &
1834	& * ( ssh (ji+jio,jj+jjo,Kmm) + gde3w(ji+jio,jj+jjo,1) ) &
1835	& / ( - gde3w(ji+jio,jj+jjo,2) + gde3w(ji+jio,jj+jjo,1) ) )
1836	END_2D
1837
1838	!--------------------------------------------------------------
1839	! 4.3 Interior faces, compute and store
1840	!-------------------------------------------------------------
1841
1842	DO_3D( 0, 0, 0, 0, 2, jpkm1 )
1843	ztmp(ji,jj,jk) = zcoef0 * ( z_rhd_pmr(ji,jj,jk,jr) + z_rhd_pmr(ji,jj,jk-1,jr) ) &
1844	& * ( - gde3w(ji+jio,jj+jjo,jk) + gde3w(ji+jio,jj+jjo,jk-1) ) &
1845	& + z_grav_10 * ( &
1846	& ( zdrhd_k_pmr(ji,jj,jk) - zdrhd_k_pmr(ji,jj,jk-1) ) &
1847	& * ( - gde3w(ji+jio,jj+jjo,jk) + gde3w(ji+jio,jj+jjo,jk-1) - z1_12 * ( zdz_k(ji+jio,jj+jjo,jk) + zdz_k (ji+jio,jj+jjo,jk-1) ) ) &
1848	& - ( zdz_k(ji+jio,jj+jjo,jk) - zdz_k(ji+jio,jj+jjo,jk-1) ) &
1849	& * ( z_rhd_pmr(ji,jj,jk,jr) - z_rhd_pmr(ji,jj,jk-1,jr) - z1_12 * ( zdrhd_k_pmr(ji,jj,jk) + zdrhd_k_pmr(ji,jj,jk-1) ) ) &
1850	& )
1851	END_3D
1852
1853	! the force on the right face could be set equal to the average of the right face for this cell and the left face for the cell to the right
1854	! this would require an lbc_lnk call
1855
1856	! lmr stands for left minus right
1857
1858	IF ( jr == 2 ) THEN
1859	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
1860	z_lmr_k(ji,jj,jk) = ztmp(ji,jj,jk) ! values on left face;
1861	END_3D
1862	ELSE
1863	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
1864	z_lmr_k(ji,jj,jk) = z_lmr_k(ji,jj,jk) - ztmp(ji,jj,jk) ! subtract the values on the right face
1865	END_3D
1866	END IF
1867
1868	IF ( ln_dbg_hpg ) CALL dbg_3dr( '4. z_lmr_k', z_lmr_k )
1869
1870	END DO ! jr
1871
1872
1873	!----------------------------------------------------------------------------------------
1874	! 5. Calculations for upper and lower faces and the vertical integration
1875	!----------------------------------------------------------------------------------------
1876
1877	z_upp_ij(:,:) = 0._wp
1878	zhpi(:,:) = 0._wp
1879	zhpj(:,:) = 0._wp
1880
1881	DO jk = 1, jpk -1
1882
1883	IF ( ln_dbg_hpg .AND. lwp ) THEN
1884	WRITE(numout,*)
1885	WRITE(numout,*) ' jk = ', jk
1886	END IF
1887
1888	!----------------------------------------------------------------------------------------
1889	! 5.1 compute and store elementary horizontal differences zfor z_rhd_pmr arrays
1890	!----------------------------------------------------------------------------------------
1891
1892	DO_2D( 0, 0, 0, 0 )
1893	zdrhd_21(ji,jj) = z_rhd_pmr(ji,jj,jk,2) - z_rhd_pmr(ji,jj,jk,1)
1894	zdrhd_32(ji,jj) = z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2)
1895	zdrhd_43(ji,jj) = z_rhd_pmr(ji,jj,jk,4) - z_rhd_pmr(ji,jj,jk,3)
1896	END_2D
1897
1898	DO_2D( 0, 0, 0, 0 )
1899	cff_31 = MAX( 2._wp * zdrhd_21(ji,jj) * zdrhd_32(ji,jj), 0._wp )
1900	z1_cff = zdrhd_21(ji,jj) + zdrhd_32(ji,jj)
1901	zz_drhd_ij(ji,jj,1) = cff_31 / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
1902
1903	cff_42 = MAX( 2._wp * zdrhd_32(ji,jj) * zdrhd_43(ji,jj), 0._wp )
1904	z1_cff = zdrhd_32(ji,jj) + zdrhd_43(ji,jj)
1905	zz_drhd_ij(ji,jj,2) = cff_42 / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
1906	END_2D
1907
1908	!----------------------------------------------------------------------------------
1909	! 5.2 apply boundary conditions at side boundaries using 5.36-5.37
1910	!----------------------------------------------------------------------------------
1911
1912
1913	! need to check this sub-section more carefully
1914
1915	zdrhd_ij(:,:,:) = zz_drhd_ij(:,:,:)
1916
1917	IF ( j_uv == 1 ) THEN
1918
1919	DO_2D( 0, 0, 0, 0)
1920	! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj)
1921	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
1922	zdrhd_ij(ji,jj,1) = aco_bc_rhd_hor * ( z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2) ) - bco_bc_rhd_hor * zz_drhd_ij(ji,jj,2)
1923	END IF
1924	! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj)
1925	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
1926	zdrhd_ij(ji,jj,2) = aco_bc_rhd_hor * ( z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2) ) - bco_bc_rhd_hor * zz_drhd_ij(ji,jj,1)
1927	END IF
1928	END_2D
1929
1930	ELSE ! j_uv == 2
1931
1932	DO_2D( 0, 0, 0, 0)
1933	! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi)
1934	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
1935	zdrhd_ij(ji,jj,1) = aco_bc_rhd_hor * ( z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2) ) - bco_bc_rhd_hor * zz_drhd_ij(ji,jj,2)
1936	END IF
1937	! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi)
1938	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
1939	zdrhd_ij(ji,jj,2) = aco_bc_rhd_hor * ( z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2) ) - bco_bc_rhd_hor * zz_drhd_ij(ji,jj,1)
1940	END IF
1941	END_2D
1942
1943	END IF ! j_uv == 2
1944
1945	IF ( ln_dbg_hpg ) THEN
1946	CALL dbg_2dr( '5.2 zdrhd_ij(:,:,1)', zdrhd_ij(:,:,1) )
1947	CALL dbg_2dr( '5.2 zdrhd_ij(:,:,2)', zdrhd_ij(:,:,2) )
1948	END IF
1949
1950	!--------------------------------------------------------------
1951	! 5.3 Calculate integrals on lower faces
1952	!-------------------------------------------------------------
1953
1954	IF ( j_uv == 1 ) THEN
1955
1956	DO_2D( 0, 0, 0, 0 )
1957	! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height)
1958	z_low_ij(ji,jj) = zcoef0 * ( z_rhd_pmr(ji,jj,jk,3) + z_rhd_pmr(ji,jj,jk,2) ) &
1959	& * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) )
1960
1961	IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN
1962	z_low_ij(ji,jj) = z_low_ij(ji,jj) - z_grav_10 * ( &
1963	& ( zdrhd_ij(ji,jj,2) - zdrhd_ij(ji,jj,1) ) &
1964	& * ( - gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_i(ji+1,jj,jk) + zdz_i(ji,jj,jk) ) ) &
1965	& - ( zdz_i(ji+1,jj,jk) - zdz_i(ji,jj,jk) ) &
1966	& * ( z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2) - z1_12 * ( zdrhd_ij(ji,jj,2) + zdrhd_ij(ji,jj,1) ) ) &
1967	& )
1968	END IF
1969	END_2D
1970
1971	IF ( ln_dbg_hpg ) CALL dbg_2dr( '5.3 z_low_ij 1', z_low_ij )
1972
1973	ELSE ! j_uv == 2
1974
1975	DO_2D( 0, 0, 0, 0 )
1976	z_low_ij(ji,jj) = zcoef0 * ( z_rhd_pmr(ji,jj,jk,3) + z_rhd_pmr(ji,jj,jk,2) ) &
1977	& * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) )
1978
1979	IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN
1980	z_low_ij(ji,jj) = z_low_ij(ji,jj) - z_grav_10 * ( &
1981	& ( zdrhd_ij(ji,jj,2) - zdrhd_ij(ji,jj,1) ) &
1982	& * ( - gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_j(ji,jj+1,jk) + zdz_j(ji,jj,jk) ) ) &
1983	& - ( zdz_j(ji,jj+1,jk) - zdz_j(ji,jj,jk) ) &
1984	& * ( z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2) - z1_12 * ( zdrhd_ij(ji,jj,2) + zdrhd_ij(ji,jj,1) ) ) &
1985	& )
1986	END IF
1987	END_2D
1988
1989	IF ( ln_dbg_hpg ) CALL dbg_2dr( '5.3 z_low_ij 2', z_low_ij )
1990
1991	END IF ! j_uv
1992
1993	!--------------------------------------------------------------
1994	! 5.4 Integrate in the vertical (including contributions from both upper and lower and side-faces)
1995	!-------------------------------------------------------------
1996	!
1997	IF ( j_uv == 1 ) THEN
1998
1999	DO_2D( 0, 0, 0, 0 )
2000	zhpi(ji,jj) = zhpi(ji,jj) + &
2001	& ( z_lmr_k(ji,jj,jk) - ( z_low_ij(ji,jj) - z_upp_ij(ji,jj) ) ) * r1_e1u(ji,jj)
2002	! add to the general momentum trend
2003	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj)
2004	END_2D
2005
2006	IF ( ln_dbg_hpg ) CALL dbg_2dr( '5.4 zhpi', zhpi )
2007
2008	ELSE ! j_uv == 2
2009
2010	DO_2D( 0, 0, 0, 0 )
2011	zhpj(ji,jj) = zhpj(ji,jj) + &
2012	& ( z_lmr_k(ji,jj,jk) - ( z_low_ij(ji,jj) - z_upp_ij(ji,jj) ) ) * r1_e2v(ji,jj)
2013	! add to the general momentum trend
2014	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj)
2015	END_2D
2016
2017	IF ( ln_dbg_hpg ) CALL dbg_2dr( '5.4 zhpj', zhpj )
2018
2019	END IF ! j_uv
2020
2021	DO_2D( 0, 0, 0, 0 )
2022	z_upp_ij(ji,jj) = z_low_ij(ji,jj)
2023	END_2D
2024
2025	END DO ! k
2026
2027	END DO ! j_uv
2028
2029	END SUBROUTINE hpg_djr
2030
2031	!-----------------------------------------------------------------------------------------------------------------
2032
2033	SUBROUTINE ref_to_tgt_cub ( ki_off_tgt, kj_off_tgt, p_dep_tgt, p_fld_tgt, p_z_ref, p_fld_ref, kk_bot_ref, p_fld_tgt_ref)
2034
2035	INTEGER, INTENT(in) :: ki_off_tgt ! offset of points in target array in i-direction
2036	INTEGER, INTENT(in) :: kj_off_tgt ! offset of points in target array in j-direction
2037	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_dep_tgt ! depths of target profiles
2038	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_fld_tgt ! field values on the target grid
2039	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_z_ref ! heights of reference profiles
2040	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_fld_ref ! field values to be interpolated (in the vertical) on reference grid
2041	INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: kk_bot_ref ! bottom levels in the reference profile
2042	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(OUT) :: p_fld_tgt_ref ! target minus reference on target grid
2043
2044	INTEGER, DIMENSION(A2D(nn_hls),jpk) :: jk_ref_for_tgt ! reference index for interpolation to target grid; target lies between jk_ref and jk_ref-1.
2045
2046	!---------------------------------------------------------------------------------------------------------------
2047
2048	INTEGER :: ji, jj, jk ! loop indices
2049	INTEGER :: jkr
2050	REAL(wp) :: z_r_6, z_r_24 ! constants
2051
2052	REAL(wp) :: zf_a, zf_b, zf_c, zf_d
2053	REAL(wp) :: zz_ref_jkr, zz_ref_jkrm1, zeta, zetasq
2054	REAL(wp) :: zf_0, zf_1, zf_2, zf_3
2055	REAL(wp) :: zave_bc, zave_ad, zdif_cb, zdif_da
2056
2057	REAL(wp) :: zz_tgt_lcl, zfld_ref_on_tgt
2058
2059	!-----------------------------------------------------------------------------------------------------------------
2060
2061	z_r_6 = 1._wp / 6._wp
2062	z_r_24 = 1._wp / 24._wp
2063
2064	! find jk_ref_for_tgt (bounding levels on reference grid for each target point
2065	CALL loc_ref_tgt ( ki_off_tgt, kj_off_tgt, p_dep_tgt, p_z_ref, kk_bot_ref, jk_ref_for_tgt )
2066
2067	DO_3D( 0, 0, 0, 0, 1, jpk-1 )
2068	zz_tgt_lcl = - p_dep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk )
2069
2070	! it would probably be better computationally for fld_ref to have the jk index first.
2071
2072	!!! jkr >= 2 and p_fld_ref has jk = 0 available
2073
2074	jkr = jk_ref_for_tgt(ji,jj,jk)
2075	zf_a = p_fld_ref(ji,jj,jkr-2)
2076	zf_b = p_fld_ref(ji,jj,jkr-1)
2077	zf_c = p_fld_ref(ji,jj,jkr )
2078	zf_d = p_fld_ref(ji,jj,jkr+1)
2079
2080	zz_ref_jkrm1 = p_z_ref( ji, jj, jkr - 1 )
2081	zz_ref_jkr = p_z_ref( ji, jj, jkr )
2082	zeta = ( zz_tgt_lcl - 0.5_wp*(zz_ref_jkr+zz_ref_jkrm1) ) / ( zz_ref_jkr - zz_ref_jkrm1 )
2083	zetasq = zeta*zeta
2084
2085	zave_bc = 0.5_wp*(zf_b+zf_c)
2086	zave_ad = 0.5_wp*(zf_a+zf_d)
2087	zdif_cb = zf_c - zf_b
2088	zdif_da = zf_d - zf_a
2089
2090	zf_0 = 1.125_wpzave_bc - 0.125_wpzave_ad
2091	zf_1 = 1.125_wpzdif_cb - z_r_24zdif_da
2092	zf_2 = 0.5_wp*(zave_ad - zave_bc) ! corrected 12/09/2021
2093	zf_3 = z_r_6 * zdif_da - 0.5_wp*zdif_cb ! corrected 12/09/2021
2094
2095	zfld_ref_on_tgt = zf_0 + zetazf_1 + zetasq(zf_2 + zeta*zf_3)
2096
2097	! when zfld_ref_on_tgt is commented out in the next line, the results for hpg_djr should agree with those for hpg_djc.
2098	p_fld_tgt_ref(ji, jj, jk) = p_fld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) - zfld_ref_on_tgt
2099
2100	END_3D
2101
2102	RETURN
2103
2104	END SUBROUTINE ref_to_tgt_cub
2105
2106	!-----------------------------------------------------------------------------------------------------------------
2107
2108	SUBROUTINE ref_to_tgt_ccs ( ki_off_tgt, kj_off_tgt, pdep_tgt, pfld_tgt, pz_ref, pfld_ref, pdfld_k_ref, kk_bot_ref, pfld_tgt_ref)
2109
2110	INTEGER, INTENT(in) :: ki_off_tgt ! offset of points in target array in i-direction
2111	INTEGER, INTENT(in) :: kj_off_tgt ! offset of points in target array in j-direction
2112	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pdep_tgt ! depths of target profiles
2113	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pfld_tgt ! field values on the target grid
2114	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pz_ref ! heights of reference profiles
2115	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pfld_ref ! reference field values
2116	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pdfld_k_ref ! harmonic averages of vertical differences of reference field
2117	INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: kk_bot_ref ! bottom levels in the reference profile
2118	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(OUT) :: pfld_tgt_ref ! target minus reference on target grid
2119
2120	!-----------------------------------------------------------------------------------------------------------------
2121	INTEGER, DIMENSION(A2D(nn_hls),jpk) :: jk_ref_for_tgt ! reference index for interpolation to target grid
2122
2123	INTEGER :: ji, jj, jk ! DO loop indices
2124	INTEGER :: jkr
2125	REAL(wp) :: z_r_6 ! constant
2126	REAL(wp) :: zz_tgt_lcl, zz_ref_jkrm1, zz_ref_jkr, zeta, zetasq
2127	REAL(wp) :: zd_dif, zd_ave, zf_dif, zf_ave
2128	REAL(wp) :: zcoef_0, zcoef_1, zcoef_2, zcoef_3
2129	REAL(wp) :: zfld_ref_on_tgt
2130	!-----------------------------------------------------------------------------------------------------------------
2131
2132	z_r_6 = 1._wp / 6._wp
2133
2134	! find jk_ref_for_tgt (bounding levels on reference grid for each target point
2135	CALL loc_ref_tgt ( ki_off_tgt, kj_off_tgt, pdep_tgt, pz_ref, kk_bot_ref, jk_ref_for_tgt )
2136
2137	DO_3D( 0, 0, 0, 0, 1, jpk-1 )
2138	zz_tgt_lcl = - pdep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk )
2139	jkr = jk_ref_for_tgt( ji, jj, jk )
2140	zz_ref_jkrm1 = pz_ref( ji, jj, jkr - 1 )
2141	zz_ref_jkr = pz_ref( ji, jj, jkr )
2142	zeta = ( zz_tgt_lcl - 0.5_wp*(zz_ref_jkr+zz_ref_jkrm1) ) / ( zz_ref_jkr - zz_ref_jkrm1 )
2143	zetasq = zeta*zeta
2144
2145	zd_dif = pdfld_k_ref(ji,jj,jkr) - pdfld_k_ref(ji,jj,jkr-1)
2146	zd_ave = 0.5_wp * ( pdfld_k_ref(ji,jj,jkr) + pdfld_k_ref(ji,jj,jkr-1) )
2147
2148	zf_dif = pfld_ref(ji,jj,jkr) - pfld_ref(ji,jj,jkr-1)
2149	zf_ave = 0.5_wp * ( pfld_ref(ji,jj,jkr) + pfld_ref(ji,jj,jkr-1) )
2150
2151	zcoef_0 = zf_ave - 0.125_wp * zd_dif
2152	zcoef_1 = 1.5_wp * zf_dif - 0.5_wp * zd_ave
2153	zcoef_2 = 0.5_wp * zd_dif
2154	zcoef_3 = 2.0_wp * zd_ave - 2._wp * zf_dif
2155
2156	zfld_ref_on_tgt = zcoef_0 + zetazcoef_1 + zetasq ( zcoef_2 + zeta*zcoef_3 )
2157
2158	! when zfld_ref_on_tgt is commented out in the next line, the results for hpg_djr should agree with those for hpg_djc.
2159
2160	pfld_tgt_ref(ji, jj, jk) = pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) - zfld_ref_on_tgt
2161
2162	! IF ( ln_dbg_hpg .AND. lwp .AND. ji == ki_dbg_cen .AND. jj == kj_dbg_cen .AND. jk == kk_dbg_cen ) THEN
2163	! WRITE(numout,*) ' zeta, zd_dif, zd_ave, zf_dif, zf_ave = ', zeta, zd_dif, zd_ave, zf_dif, zf_ave
2164	! WRITE(numout,*) ' zz_tgt_lcl, jkr, zz_ref_jkr, zz_ref_jkrm1 =', zz_tgt_lcl, jkr, zz_ref_jkr, zz_ref_jkrm1
2165	! WRITE(numout,*) ' pfld_ref(ji,jj,jkr), pfld_ref(ji,jj,jkr-1) = ', pfld_ref(ji,jj,jkr), pfld_ref(ji,jj,jkr-1)
2166	! WRITE(numout,*) ' zfld_ref_on_tgt = ', zfld_ref_on_tgt
2167	! WRITE(numout,*) ' pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) = ', pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk)
2168	! END IF
2169
2170	END_3D
2171
2172	IF ( ln_dbg_hpg ) CALL dbg_3dr( 'ref_to_tgt_ccs: pfld_tgt_ref', pfld_tgt_ref )
2173
2174	RETURN
2175
2176	END SUBROUTINE ref_to_tgt_ccs
2177
2178	!-----------------------------------------------------------------------------------------------------------------
2179
2180	SUBROUTINE loc_ref_tgt ( ki_off_tgt, kj_off_tgt, p_dep_tgt, p_z_ref, kk_bot_ref, kk_ref_for_tgt )
2181
2182	INTEGER, INTENT(in) :: ki_off_tgt ! offset of points in target array in i-direction
2183	INTEGER, INTENT(in) :: kj_off_tgt ! offset of points in target array in j-direction
2184	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_dep_tgt ! depths of target profiles
2185	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_z_ref ! heights of reference profiles
2186	INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: kk_bot_ref ! bottom levels in the reference profile
2187	INTEGER, DIMENSION(A2D(nn_hls),jpk), INTENT(OUT) :: kk_ref_for_tgt ! reference index for interpolation to target grid (the lower point)
2188	! with jkr = kk_ref_for_tgt(ji,jj,jk) the reference levels are jkr-1 and jkr
2189	!-----------------------------------------------------------------------------------------------------------------
2190
2191	INTEGER, DIMENSION(A2D(nn_hls)) :: jk_tgt, jk_ref ! vertical level being processed on target and reference grids respectively
2192
2193	INTEGER :: ji, jj, jk_comb, jk_bot
2194
2195	INTEGER :: jk, jkr, jcount ! for debugging only
2196	REAL(wp):: z_tgt, z_below, z_above ! for debugging only
2197
2198	INTEGER :: jk_init ! initial jk value for search
2199
2200	REAL(wp):: tol_wp ! this should be the working precision tolerance
2201
2202	tol_wp = 1.E-4 ! the inferred bottom depth seems to be accurate to about 1.E-6.
2203
2204	!-----------------------------------------------------------------------------------------------------------------
2205
2206	! 1. Initialisation for search for points on reference grid bounding points on the target grid
2207	! the first point on the target grid is assumed to lie between the first two points on the reference grid
2208	IF ( ln_hpg_djr_ref_ccs ) THEN
2209	jk_init = 2
2210	ELSE
2211	jk_init = 3 ! for simple cubic use off-centred interpolation near the surface
2212	END IF
2213
2214	jk_ref(:,:) = jk_init
2215	kk_ref_for_tgt(:,:,1) = jk_init
2216	jk_tgt(:,:) = 2
2217
2218	! 2. Find kk_ref_for_tgt
2219
2220	DO jk_comb = 1, 2*(jpk-1)
2221
2222	! if level jk_tgt in target profile is lower than jk_ref in reference profile add one to jk_ref
2223	DO_2D( 0, 0, 0, 0 )
2224	IF ( - p_dep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk_tgt(ji,jj) ) < p_z_ref( ji, jj, jk_ref(ji,jj) ) ) THEN
2225	IF ( jk_ref(ji,jj) < jpk-1 ) jk_ref(ji,jj) = jk_ref(ji,jj) + 1
2226	END IF
2227	END_2D
2228
2229	! if level jk_tgt lies above or at same level as jk_ref, store jk_ref for this level and add one to jk_tgt
2230	DO_2D( 0, 0, 0, 0 )
2231	IF ( - p_dep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk_tgt(ji,jj) ) + tol_wp > p_z_ref( ji, jj, jk_ref(ji,jj) ) ) THEN
2232	IF ( jk_tgt(ji,jj) < jpk ) THEN
2233	kk_ref_for_tgt( ji, jj, jk_tgt(ji,jj) ) = jk_ref(ji,jj)
2234	jk_tgt(ji,jj) = jk_tgt(ji,jj) + 1
2235	END IF
2236	END IF
2237	END_2D
2238
2239	IF ( lwp .AND. ln_dbg_hpg ) THEN
2240	CALL dbg_2di_k( 'jk_ref', jk_ref, jk_comb )
2241	CALL dbg_2di_k( 'jk_tgt', jk_tgt, jk_comb )
2242	END IF
2243
2244	END DO ! jk_comb
2245
2246
2247	! 3. Checks to make sure we do not read or write out of bounds
2248
2249	! 3.1 First test on jk_tgt to check that it reaches the bottom level mbkt
2250
2251	jcount = 0
2252	DO_2D(0,0,0,0)
2253	IF ( jk_tgt(ji,jj) < (mbkt(ji+ki_off_tgt, jj+kj_off_tgt) + 1) ) jcount = jcount + 1
2254	END_2D
2255
2256	IF ( jcount > 0 ) THEN
2257	WRITE( numout,*) 'loc_ref_tgt: stopping because kk_ref_for_tgt will not cover all sea-points; jcount = ', jcount
2258	CALL ctl_stop( 'dyn_hpg_djr : kk_ref_for_tgt will not cover all sea-points' )
2259	END IF
2260
2261	! 3.2 kk_ref_for_tgt pointing to invalid level in reference profile
2262
2263	jcount = 0
2264	DO_2D(0,0,0,0)
2265	jk_bot = mbkt(ji+ki_off_tgt, jj+kj_off_tgt)
2266	IF ( kk_ref_for_tgt (ji,jj,jk_bot) > kk_bot_ref(ji,jj) ) jcount = jcount + 1
2267	END_2D
2268
2269	IF ( jcount > 0 ) THEN
2270	WRITE( numout,*) 'loc_ref_tgt: stopping because kk_ref_for_tgt points to a level below the bottom of the reference profile; jcount = ', jcount
2271	CALL ctl_stop( 'dyn_hpg_djr : kk_ref_for_tgt points to a level below the bottom of the reference profile' )
2272	END IF
2273
2274
2275	! 3.3 diagnostic check that the right reference levels are chosen for the target profile
2276	IF ( lwp .AND. ln_dbg_hpg ) THEN
2277	WRITE( numout,*) 'loc_ref_tgt: ki_off_tgt, kj_off_tgt = ', ki_off_tgt, kj_off_tgt
2278	CALL dbg_3dr( '-p_dep_tgt', -p_dep_tgt )
2279	CALL dbg_3dr( 'p_z_ref', p_z_ref )
2280	CALL dbg_3di( 'kk_ref_for_tgt', kk_ref_for_tgt )
2281
2282	jcount = 0
2283	DO_3D(0,0,0,0,2,jpk-1)
2284	z_tgt = - p_dep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk)
2285	jkr = kk_ref_for_tgt(ji,jj,jk)
2286	z_below = p_z_ref( ji, jj, jkr )
2287	z_above = p_z_ref( ji, jj, jkr-1 )
2288	IF ( ( z_above < z_tgt .AND. jkr > jk_init ) .OR. z_below > (z_tgt + tol_wp) ) THEN
2289	IF ( jcount < 10 ) THEN
2290	WRITE(numout,*) 'loc_ref_tgt: ji, jj, jk, ki_off_tgt, kj_off_tgt, z_tgt, z_below, z_above '
2291	WRITE(numout,*) ji, jj, jk, ki_off_tgt, kj_off_tgt, z_tgt, z_below, z_above
2292	END IF
2293	jcount = jcount + 1
2294	END IF
2295	END_3D
2296	WRITE(numout,*) 'loc_ref_tgt: jcount = ', jcount
2297	END IF
2298
2299	! 3.4 Same check as in 4.2 but detailed diagnostics not written out.
2300	jcount = 0
2301	DO_3D(0,0,0,0,2,jpk-1)
2302	z_tgt = - p_dep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk)
2303	jkr = kk_ref_for_tgt(ji,jj,jk)
2304	z_below = p_z_ref( ji, jj, jkr )
2305	z_above = p_z_ref( ji, jj, jkr-1 )
2306	IF ( ( z_above < z_tgt .AND. jkr > jk_init) .OR. z_below > (z_tgt + tol_wp) ) THEN
2307	jcount = jcount + 1
2308	END IF
2309	END_3D
2310
2311	IF ( jcount > 0 ) THEN
2312	WRITE( numout,*) 'loc_ref_tgt: stopping because jcount is non-zero; jcount = ', jcount
2313	CALL ctl_stop( 'dyn_hpg_djr : loc_ref_tgt failed to locate target points correctly' )
2314	END IF
2315
2316	! 4. Adjust kk_ref_for_tgt so that interpolation of simple cubic is off-centred at the bottom (does not require boundary conditions)
2317	! This assumes that every sea point has at least 4 levels.
2318
2319	IF ( .NOT. ln_hpg_djr_ref_ccs ) THEN
2320	DO_3D(0, 0, 0, 0, 2, jpk-1)
2321	IF ( kk_ref_for_tgt(ji,jj, jk) == kk_bot_ref(ji,jj) ) kk_ref_for_tgt(ji,jj, jk) = kk_ref_for_tgt(ji,jj, jk) - 1
2322	END_3D
2323	END IF
2324
2325	RETURN
2326
2327	END SUBROUTINE loc_ref_tgt
2328
2329	!----------------------------------------------------------------------------
2330
2331	SUBROUTINE hpg_ffr( kt, Kmm, puu, pvv, Krhs )
2332	!!---------------------------------------------------------------------
2333	!! * ROUTINE hpg_ffr *
2334	!!
2335	!! ** Method : s-coordinate case forces on faces scheme.
2336	!! Local density subtracted using cubic or constrained cubic splines (ccs)
2337	!! Remaining densities interpolated to centre either linearly or with ccs
2338	!! Vertical representation either using quadratic density or classic 2nd order accurate
2339	!!
2340	!! ** Action : - Update puu(..,Krhs) and pvv(..,Krhs) with the now hydrastatic pressure trend
2341	!!----------------------------------------------------------------------
2342	INTEGER , INTENT( in ) :: kt ! ocean time-step index
2343	INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices
2344	REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
2345	!!
2346
2347	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zrhd_ref, zdrhd_k_ref ! densities (rhd) of reference profile
2348	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zz_ref ! heights of reference profile
2349	REAL(wp), DIMENSION(A2D(nn_hls),jpk,4) :: z_rhd_pmr ! profiles minus reference
2350
2351	! The following fields could probably be made single level or at most 2 level fields
2352	REAL(wp), DIMENSION(A2D(nn_hls),jpk,4) :: z_e3t_pmr, z_depw_pmr ! corresponding e3t and gdepw pmr profiles
2353	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z_rhd_mid, z_e3t_mid ! profiles and e3t interpolated to middle of cell (using ccs)
2354	REAL(wp), DIMENSION(A2D(nn_hls),jpk,3) :: z_ddepw_ij ! constrained spline horizontal differences in gdepw
2355	REAL(wp), DIMENSION(A2D(nn_hls),jpk,3) :: z_p_pmr, z_F_pmr ! pressures and forces on vertical faces
2356	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z_F_upp ! forces on upper faces
2357
2358	INTEGER, DIMENSION(A2D(nn_hls)) :: jk_bot_ref ! bottom level in reference profiles
2359	INTEGER, DIMENSION(A2D(nn_hls),3) :: j_mbkt ! bottom levels for the 3 profiles in the cell
2360
2361	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zpforces ! total of forces
2362
2363	INTEGER :: ji, jj, jk, j_uv, jr ! dummy loop indices
2364	INTEGER :: jio, jjo ! offsets for the 2 point stencil (0 or 1)
2365	INTEGER :: ji_ro, jj_ro, jr_offset ! offsets for the reference profile
2366	INTEGER :: jn_hor_pts ! number of profiles required in horizontal (4 if cubic interpolarion or 2 if not)
2367
2368	REAL(wp) :: z_r_6 ! 1/6
2369	REAL(wp) :: zr_a, zr_b, zr_c ! rhd evaluated at i, i+1/2 and i+1
2370	REAL(wp) :: za_0, za_1, za_2 ! polynomial coefficients
2371	!!----------------------------------------------------------------------
2372
2373	IF( kt == nit000 ) THEN
2374	IF(lwp) WRITE(numout,*)
2375	IF(lwp) WRITE(numout,*) 'dyn:hpg_Lin : hydrostatic pressure gradient trend'
2376	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, forces on faces with reference removed'
2377	ENDIF
2378	!
2379
2380	z_r_6 = 1.0_wp/6.0_wp
2381
2382	IF ( ln_hpg_ffr_hor_ccs .OR. ln_hpg_ffr_hor_cub ) THEN
2383	jn_hor_pts = 4
2384	jr_offset = 2
2385	ELSE
2386	jn_hor_pts = 2
2387	jr_offset = 1
2388	END IF
2389
2390	! WRITE( numout, *) ' hpg_ffr: kt, jn_hor_pts, jr_offset = ', kt, jn_hor_pts, jr_offset
2391
2392	DO j_uv = 1, 2
2393
2394	IF ( j_uv == 1 ) THEN
2395	jio = 1 ! set for the whole of the j_uv loop (loop ends near the bottom of the routine)
2396	jjo = 0
2397	ELSE
2398	jio = 0
2399	jjo = 1
2400	END IF
2401
2402	! 1. Calculate the referene profile from all points in the stencil
2403
2404	IF ( ln_hpg_ffr_ref ) THEN
2405	CALL calc_rhd_ref(j_uv, jn_hor_pts, zrhd_ref, zz_ref, jk_bot_ref) ! Uses rhd (IN) to calculate all other fields (OUT)
2406	IF ( ln_hpg_ffr_ref_ccs ) CALL calc_drhd_k(zrhd_ref, jk_bot_ref, zdrhd_k_ref)
2407	END IF
2408
2409	IF ( ln_dbg_hpg ) THEN
2410	CALL dbg_3dr( '2. rhd', rhd )
2411	CALL dbg_3dr( '2. gdept', gdept(:,:,:,Kmm) )
2412	END IF
2413
2414	! 2. Interpolate reference profile to all points in the stencil and calculate z_rhd_pmr (profile minus reference)
2415
2416	DO jr = 1, jn_hor_pts
2417
2418
2419	IF ( j_uv == 1 ) THEN
2420	ji_ro = jr - jr_offset ! range of jio: is -1 to 2 for jn_hor_pts = 4; is 0 to 1 for jn_hor_pts = 2
2421	jj_ro = 0
2422	ELSE
2423	ji_ro = 0
2424	jj_ro = jr - jr_offset
2425	END IF
2426
2427	IF ( ln_hpg_ffr_ref ) THEN !!!! MJB should gde3w below now be gdept ?? !!!!
2428
2429	IF ( ln_hpg_ffr_ref_ccs ) THEN
2430	CALL ref_to_tgt_ccs ( ji_ro, jj_ro, gde3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )
2431	ELSE
2432	CALL ref_to_tgt_cub ( ji_ro, jj_ro, gde3w, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )
2433	END IF
2434
2435	ELSE ! no subtraction of reference profile
2436
2437	DO_3D (0,0,0,0,1,jpk-1)
2438	z_rhd_pmr(ji,jj,jk,jr) = rhd(ji+ji_ro,jj+jj_ro,jk)
2439	END_3D
2440
2441	END IF ! ln_hpg_ffr_ref ) THEN
2442
2443	DO_3D (0,0,0,0,1,jpk)
2444	z_e3t_pmr (ji,jj,jk,jr) = e3t(ji+ji_ro,jj+jj_ro,jk,Kmm)
2445	z_depw_pmr(ji,jj,jk,jr) = gdepw(ji+ji_ro,jj+jj_ro,jk,Kmm)
2446	END_3D
2447
2448	IF ( ln_dbg_hpg ) THEN
2449	CALL dbg_3dr( '2. z_rhd_pmr', z_rhd_pmr(:,:,:,jr) )
2450	END IF
2451
2452	END DO
2453
2454	! 3. Do horizontal interpolation to form intermediate densities (either linear or cubic or constrained cubic)
2455	! Transfers data points from reference stencil (2 or 4 point) to a 3 point horizontal stencil
2456
2457	IF ( ln_hpg_ffr_hor_ccs .OR. ln_hpg_ffr_hor_cub) THEN
2458
2459	IF ( ln_hpg_ffr_hor_ccs ) THEN
2460	CALL calc_mid_ccs(j_uv, aco_bc_rhd_hor, bco_bc_rhd_hor, z_rhd_pmr, z_rhd_mid)
2461	CALL calc_mid_ccs(j_uv, aco_bc_z_hor, bco_bc_z_hor, z_e3t_pmr, z_e3t_mid)
2462	CALL calc_dz_dij_ccs( j_uv, z_depw_pmr, z_ddepw_ij )
2463	ELSE ! ln_hpg_ffr_hor_cub
2464	CALL calc_mid_cub(j_uv, aco_bc_rhd_hor, bco_bc_rhd_hor, z_rhd_pmr, z_rhd_mid)
2465	CALL calc_mid_cub(j_uv, aco_bc_z_hor, bco_bc_z_hor, z_e3t_pmr, z_e3t_mid)
2466	CALL calc_dz_dij_cub( j_uv, z_depw_pmr, z_ddepw_ij )
2467	END IF
2468
2469	DO_3D (0,0,0,0,1,jpk-1)
2470	z_rhd_pmr(ji,jj,jk,1) = z_rhd_pmr(ji,jj,jk,2)
2471	z_rhd_pmr(ji,jj,jk,2) = z_rhd_mid(ji,jj,jk)
2472	z_e3t_pmr(ji,jj,jk,1) = z_e3t_pmr(ji,jj,jk,2)
2473	z_e3t_pmr(ji,jj,jk,2) = z_e3t_mid(ji,jj,jk)
2474	END_3D
2475
2476	ELSE ! simple linear interpolation
2477
2478	DO_3D (0,0,0,0,1,jpk-1)
2479	z_rhd_pmr(ji,jj,jk,3) = z_rhd_pmr(ji,jj,jk,2)
2480	z_rhd_pmr(ji,jj,jk,2) = 0.5_wp*( z_rhd_pmr(ji,jj,jk,1) + z_rhd_pmr(ji,jj,jk,2) )
2481	END_3D
2482	DO_3D (0,0,0,0,1,jpk-1)
2483	z_e3t_pmr(ji,jj,jk,1) = e3t(ji, jj, jk, Kmm)
2484	z_e3t_pmr(ji,jj,jk,2) = 0.5_wp*( e3t(ji, jj, jk, Kmm) + e3t(ji+jio, jj+jjo, jk, Kmm) )
2485	z_e3t_pmr(ji,jj,jk,3) = e3t(ji+jio, jj+jjo, jk, Kmm)
2486	END_3D
2487	END IF
2488
2489	DO_2D (0,0,0,0)
2490	j_mbkt(ji,jj,1) = mbkt(ji, jj)
2491	j_mbkt(ji,jj,2) = MIN( mbkt(ji, jj), mbkt(ji+jio, jj+jjo) )
2492	j_mbkt(ji,jj,3) = mbkt(ji+jio, jj+jjo)
2493	END_2D
2494
2495	IF ( ln_dbg_hpg ) THEN
2496	CALL dbg_3dr( '3. e3t ', e3t )
2497	CALL dbg_3dr( '3. z_rhd_pmr: 1', z_rhd_pmr(:,:,:,1) )
2498	CALL dbg_3dr( '3. z_rhd_pmr: 2', z_rhd_pmr(:,:,:,2) )
2499	CALL dbg_3dr( '3. z_rhd_pmr: 3', z_rhd_pmr(:,:,:,3) )
2500	END IF
2501
2502	! 4. vertical interpolation of densities, calculating pressures and forces on vertical faces between w levels
2503
2504	IF ( ln_hpg_ffr_vrt_quad ) THEN
2505	DO jr = 1, 3
2506	CALL vrt_int_quad(j_mbkt(:,:,jr), z_rhd_pmr(:,:,:,jr), z_e3t_pmr(:,:,:,jr), z_p_pmr(:,:,:,jr), z_F_pmr(:,:,:,jr))
2507	END DO ! jr
2508
2509	ELSE ! .NOT. ln_hpg_ffr_vrt_quad (simplest vertical integration scheme)
2510
2511	DO jr = 1, 3
2512	DO_2D(0,0,0,0)
2513	z_p_pmr(ji,jj,1,jr) = 0._wp
2514	END_2D
2515	DO jk = 1, jpk - 1
2516	DO_2D(0,0,0,0)
2517	z_p_pmr(ji,jj,jk+1,jr) = z_p_pmr(ji,jj,jk,jr) + gravz_e3t_pmr(ji,jj,jk,jr)z_rhd_pmr(ji,jj,jk,jr)
2518	END_2D
2519	DO_2D(0,0,0,0)
2520	z_F_pmr(ji,jj,jk,jr) = 0.5_wpz_e3t_pmr(ji,jj,jk,jr)( z_p_pmr(ji,jj,jk,jr) + z_p_pmr(ji,jj,jk+1,jr) )
2521	END_2D
2522	END DO ! jk
2523	END DO ! jr
2524
2525	END IF ! ln_hpg_ffr_vrt_quad
2526
2527	IF ( ln_dbg_hpg ) THEN
2528	CALL dbg_3dr( '4. z_p_pmr: 1', z_p_pmr(:,:,:,1) )
2529	CALL dbg_3dr( '4. z_p_pmr: 2', z_p_pmr(:,:,:,2) )
2530	CALL dbg_3dr( '4. z_p_pmr: 3', z_p_pmr(:,:,:,3) )
2531	CALL dbg_3dr( '4. z_F_pmr: 1', z_F_pmr(:,:,:,1) )
2532	CALL dbg_3dr( '4. z_F_pmr: 3', z_F_pmr(:,:,:,3) )
2533	CALL dbg_3dr( ' z_e3t_pmr: 1', z_e3t_pmr(:,:,:,1) )
2534	CALL dbg_3dr( ' z_e3t_pmr: 2', z_e3t_pmr(:,:,:,2) )
2535	CALL dbg_3dr( ' z_e3t_pmr: 3', z_e3t_pmr(:,:,:,3) )
2536	END IF
2537
2538
2539	! 5. Calculate forces on the upper faces and hence on the total forces on the cells (zpforces)
2540
2541	DO_2D(0,0,0,0)
2542	z_F_upp(ji,jj,1) = 0._wp
2543	END_2D
2544
2545	IF ( ln_hpg_ffr_hor_ccs .OR. ln_hpg_ffr_hor_cub ) THEN ! use Simpson's rule
2546	DO_3D(0,0,0,0,2,jpk)
2547	z_F_upp(ji,jj,jk) = - z_r_6 * ( z_ddepw_ij(ji,jj,jk,1)*z_p_pmr(ji,jj,jk,1) &
2548	& + 4._wpz_ddepw_ij(ji,jj,jk,2)z_p_pmr(ji,jj,jk,2) &
2549	& + z_ddepw_ij(ji,jj,jk,3)*z_p_pmr(ji,jj,jk,3) )
2550	END_3D
2551	ELSE ! use trapezoidal rule
2552	DO_3D(0,0,0,0,2,jpk)
2553	z_F_upp(ji,jj,jk) = 0.5_wp * ( gdepw(ji,jj,jk,Kmm) - gdepw(ji+jio,jj+jjo,jk,Kmm) ) * ( z_p_pmr(ji,jj,jk,1) + z_p_pmr(ji,jj,jk,3) )
2554	&
2555	END_3D
2556	END IF
2557
2558	IF ( ln_dbg_hpg ) THEN
2559	CALL dbg_3dr( '4. z_F_upp: ', z_F_upp )
2560	CALL dbg_3dr( '4. gdepw: ', gdepw )
2561	CALL dbg_3dr( '4. z_ddepw_ij 1: ', z_ddepw_ij(:,:,:,1) )
2562	CALL dbg_3dr( '4. z_ddepw_ij 2: ', z_ddepw_ij(:,:,:,2) )
2563	CALL dbg_3dr( '4. z_ddepw_ij 3: ', z_ddepw_ij(:,:,:,3) )
2564
2565	END IF
2566
2567	IF ( j_uv == 1 ) THEN
2568	DO_3D(0,0,0,0,1,jpk-1)
2569	zpforces(ji,jj,jk) = z_F_pmr(ji,jj,jk,1) - z_F_pmr(ji,jj,jk,3) + z_F_upp(ji,jj,jk) - z_F_upp(ji,jj,jk+1)
2570	puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zpforces(ji,jj,jk) / (e1u(ji,jj)*e3u(ji,jj,jk,Kmm))
2571	END_3D
2572	IF ( ln_dbg_hpg ) CALL dbg_3dr( '5. u zpforces: ', zpforces )
2573
2574	ELSE
2575	DO_3D(0,0,0,0,1,jpk-1)
2576	zpforces(ji,jj,jk) = z_F_pmr(ji,jj,jk,1) - z_F_pmr(ji,jj,jk,3) + z_F_upp(ji,jj,jk) - z_F_upp(ji,jj,jk+1)
2577	pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zpforces(ji,jj,jk) / (e1v(ji,jj)*e3v(ji,jj,jk,Kmm))
2578	END_3D
2579	IF ( ln_dbg_hpg ) CALL dbg_3dr( '5. v zpforces: ', zpforces )
2580	END IF
2581
2582	! temporary output of fields for debugging etc.
2583	IF ( j_uv == 1) THEN
2584	CALL iom_put( "gdepw_hpg", gdepw(:,:,:,Kmm) )
2585	CALL iom_put( "rhd_hpg", rhd )
2586	CALL iom_put( "pressure", z_p_pmr(:,:,:,1) )
2587	CALL iom_put( "u_force_west", z_F_pmr(:,:,:,1) )
2588	CALL iom_put( "u_force_upper", z_F_upp )
2589	ELSE
2590	CALL iom_put( "v_force_south", z_F_pmr(:,:,:,1) )
2591	CALL iom_put( "v_force_upper", z_F_upp )
2592	END IF
2593
2594	END DO ! j_uv
2595	!
2596	END SUBROUTINE hpg_ffr
2597
2598	!------------------------------------------------------------------------------------------
2599
2600	SUBROUTINE calc_rhd_ref (k_uv, kn_hor, prhd_ref, pz_ref, kk_bot_ref)
2601
2602	!!---------------------------------------------------------------------
2603	!! * ROUTINE calc_rhd_ref *
2604	!!
2605	!! ** Method : Find the deepest cell within the stencil.
2606	!! (Later will extend to producing a reference profile that spans the highest and lowest points in the stencil)
2607	!!
2608	!!
2609	!!
2610	!! ** Action : - Set prhd_ref, pz_ref, kk_bot_ref
2611	!!----------------------------------------------------------------------
2612	INTEGER , INTENT( in ) :: k_uv ! 1 for u-vel; 2 for v_vel
2613	INTEGER , INTENT( in ) :: kn_hor ! 4 for cubic; 2 for linear
2614
2615	REAL(wp), DIMENSION(:,:,:), INTENT(out) :: prhd_ref ! densities of reference profile
2616	REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pz_ref ! heights of " "
2617	INTEGER, DIMENSION(:,:), INTENT(out) :: kk_bot_ref ! bottom level of " "
2618
2619	INTEGER, DIMENSION(A2D(nn_hls)) :: ji_ref, jj_ref
2620
2621	INTEGER ji, jj, jk ! standard loop indices
2622	INTEGER jio, jjo ! offsets
2623	INTEGER jib, jjb ! second set of indices to deeper points
2624	INTEGER jir, jjr ! indices of reference profile
2625	REAL(wp) zhta, zhtb
2626
2627	IF (k_uv == 1) THEN
2628	jio = 1
2629	jjo = 0
2630	ELSE
2631	jio = 0
2632	jjo = 1
2633	END IF
2634
2635
2636	DO_2D( 0, 0, 0, 0 )
2637
2638	IF ( ht_0(ji+jio,jj+jjo) >= ht_0(ji,jj) ) THEN
2639	zhta = ht_0(ji+jio,jj+jjo)
2640	ji_ref(ji,jj) = ji+jio
2641	jj_ref(ji,jj) = jj+jjo
2642	ELSE
2643	zhta = ht_0(ji,jj)
2644	ji_ref(ji,jj) = ji
2645	jj_ref(ji,jj) = jj
2646	END IF
2647
2648	IF ( kn_hor == 4 ) THEN
2649	IF ( ht_0(ji-jio,jj-jjo) >= ht_0(ji+2jio,jj+2jjo) ) THEN
2650	zhtb = ht_0(ji-jio,jj-jjo)
2651	jib = ji-jio
2652	jjb = jj-jjo
2653	ELSE
2654	zhtb = ht_0(ji+2jio,jj+2jjo)
2655	jib = ji+2*jio
2656	jjb = jj+2*jjo
2657	END IF
2658	IF ( zhta < zhtb ) THEN
2659	ji_ref(ji,jj) = jib
2660	jj_ref(ji,jj) = jjb
2661	END IF
2662	END IF
2663
2664	END_2D
2665
2666	DO_3D( 0, 0, 0, 0, 1, jpk-1)
2667	jir = ji_ref(ji,jj)
2668	jjr = jj_ref(ji,jj)
2669	prhd_ref (ji,jj,jk) = rhd(jir, jjr, jk)
2670	pz_ref (ji,jj,jk) = - gde3w(jir, jjr, jk)
2671	END_3D
2672
2673	DO_2D( 0, 0, 0, 0)
2674	jir = ji_ref(ji,jj)
2675	jjr = jj_ref(ji,jj)
2676	kk_bot_ref(ji,jj) = mbkt(jir,jjr)
2677	END_2D
2678
2679	END SUBROUTINE calc_rhd_ref
2680
2681	!------------------------------------------------------------------------------------------
2682
2683	SUBROUTINE calc_drhd_k(p_rhd, kk_bot, p_drhd_k)
2684
2685	!!---------------------------------------------------------------------
2686	!! * ROUTINE calc_drhd_k *
2687	!!
2688	!! ** Method : Calculate harmonic averages of vertical differences and apply upper and lower boundary conditions
2689	!!
2690	!! ** Action : - Set p_drhd_k
2691	!!----------------------------------------------------------------------
2692
2693	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_rhd ! densities of profile
2694	INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: kk_bot ! bottom level of profile
2695	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(out) :: p_drhd_k ! harmonic mean of vertical differences of profile
2696
2697	INTEGER ji, jj, jk ! standard loop indices
2698	INTEGER jio, jjo ! offsets
2699	INTEGER iktb ! index of the bottom of ref profile
2700	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: ztmp ! primitive vertical differences
2701
2702	REAL(wp) cffw, z1_cff, zep
2703
2704	DO_3D( 0, 0, 0, 0, 2, jpk )
2705	ztmp(ji,jj,jk) = p_rhd(ji,jj,jk) - p_rhd(ji,jj,jk-1)
2706	END_3D
2707
2708	zep = 1.e-15
2709	DO_3D( 0, 0, 0, 0, 2, jpkm1 )
2710	cffw = MAX( 2._wp * ztmp(ji,jj,jk) * ztmp(ji,jj,jk+1), 0._wp )
2711	z1_cff = ztmp(ji,jj,jk) + ztmp(ji,jj,jk+1)
2712	p_drhd_k(ji,jj,jk) = cffw / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
2713	END_3D
2714
2715	DO_2D( 0, 0, 0, 0 )
2716	p_drhd_k(ji,jj,1) = aco_bc_rhd_srf * ( p_rhd(ji,jj,2) - p_rhd(ji,jj,1) ) - bco_bc_rhd_srf * p_drhd_k(ji,jj,2)
2717	iktb = kk_bot(ji,jj)
2718	IF ( iktb > 1 ) THEN
2719	p_drhd_k(ji,jj,iktb) = aco_bc_rhd_bot * (p_rhd(ji,jj,iktb) - p_rhd(ji,jj,iktb-1) ) - bco_bc_rhd_bot * p_drhd_k(ji,jj,iktb-1)
2720	END IF
2721	END_2D
2722
2723	END SUBROUTINE calc_drhd_k
2724
2725	!----------------------------------------------------------------------------
2726
2727	SUBROUTINE calc_mid_ccs( k_uv, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_fld_mid )
2728
2729	!!---------------------------------------------------------------------
2730	!! * ROUTINE calc_mid_ccs *
2731	!!
2732	!! ** Method : Use constrained cubic spline to interpolate 4 profiles to their central point
2733	!! This version can only be used to interpolate fields on tracer levels (not w-levels)
2734	!!
2735	!! ** Action : - set p_fld_mid
2736	!!----------------------------------------------------------------------
2737
2738	INTEGER , INTENT(in) :: k_uv ! 1 for u-vel; 2 for v_vel
2739	REAL(wp) , INTENT(in) :: p_aco_bc_fld_hor ! a coeff for horizontal bc (von Neumann or linear extrapolation)
2740	REAL(wp) , INTENT(in) :: p_bco_bc_fld_hor ! b coeff for horizontal bc (von Neumann or linear extrapolation)
2741	REAL(wp), DIMENSION(A2D(nn_hls),jpk,4), INTENT(in) :: p_fld_pmr ! field in pmr form
2742	REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(out) :: p_fld_mid ! field at mid-point
2743
2744	REAL(wp), DIMENSION(A2D(nn_hls),jpk,3) :: zz_dfld_ij ! constrained spline horizontal differences (dimension 3 for consistency with calc_dfld_cen_dij)
2745
2746	INTEGER ji, jj, jk ! standard loop indices
2747	INTEGER :: j_t_levs ! indicator that field passed is on tracer levels
2748
2749	j_t_levs = 1
2750
2751	DO jk = 1, jpk
2752
2753	CALL calc_dfld_pmr_ij( k_uv, j_t_levs, jk, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, zz_dfld_ij )
2754
2755	! first contribution is simple centred average; p_rhd_pmr has 4 points; 2 and 3 are the central ones
2756	DO_2D( 0, 0, 0, 0)
2757	p_fld_mid(ji,jj,jk) = 0.5_wp * ( p_fld_pmr(ji,jj,jk,2) + p_fld_pmr(ji,jj,jk,3) )
2758	END_2D
2759
2760	IF ( k_uv == 1 ) THEN
2761	DO_2D( 0, 0, 0, 0)
2762	IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN
2763	p_fld_mid(ji,jj,jk) = p_fld_mid(ji,jj,jk) - 0.125_wp * ( zz_dfld_ij(ji,jj,jk,2) - zz_dfld_ij(ji,jj,jk,1) )
2764	END IF
2765	END_2D
2766	ELSE ! k_uv == 2
2767	DO_2D( 0, 0, 0, 0)
2768	IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN
2769	p_fld_mid(ji,jj,jk) = p_fld_mid(ji,jj,jk) - 0.125_wp * ( zz_dfld_ij(ji,jj,jk,2) - zz_dfld_ij(ji,jj,jk,1) )
2770	END IF
2771	END_2D
2772	END IF ! k_uv
2773
2774	END DO ! jk
2775
2776	END SUBROUTINE calc_mid_ccs
2777
2778	!----------------------------------------------------------------------------
2779
2780	SUBROUTINE calc_mid_cub( k_uv, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_fld_mid )
2781
2782	!!---------------------------------------------------------------------
2783	!! * ROUTINE calc_mid_cub *
2784	!!
2785	!! ** Method : Use simple cubic polynomial to interpolate 4 profiles to their central point
2786	!! The coefficients p_aco_bc_fld_hor, p_bco_bc_fld_hor are not currently used.
2787	!! The simplest form of von Neumann conditions horizontal bc is used (one could off-centred polynomials)
2788	!! ** Action : - set p_fld_mid
2789	!!----------------------------------------------------------------------
2790
2791	INTEGER , INTENT(in) :: k_uv ! 1 for u-vel; 2 for v_vel
2792	REAL(wp) , INTENT(in) :: p_aco_bc_fld_hor ! a coeff for horizontal bc (von Neumann or linear extrapolation)
2793	REAL(wp) , INTENT(in) :: p_bco_bc_fld_hor ! b coeff for horizontal bc (von Neumann or linear extrapolation)
2794	REAL(wp), DIMENSION(A2D(nn_hls),jpk,4), INTENT(inout) :: p_fld_pmr ! field in pmr form
2795	REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(out) :: p_fld_mid ! field at mid-point
2796
2797	REAL(wp), DIMENSION(A2D(nn_hls)) :: zdfld_21, zdfld_32, zdfld_43 ! primitive horizontal differences
2798	REAL(wp), DIMENSION(A2D(nn_hls),3) :: zz_dfld_ij ! constrained spline horizontal differences (dimension 3 for consistency with calc_dfld_cen_dij)
2799
2800	INTEGER ji, jj, jk ! standard loop indices
2801	REAL(wp) z_1_16, z_9_16 ! temporary sum and products
2802	REAL(wp) z_cor ! correction to the central value
2803
2804	z_1_16 = 1._wp / 16._wp ; z_9_16 = 9._wp / 16._wp
2805
2806	DO jk = 1, jpk
2807
2808	! first contribution is simple centred average plus 1/16 of it; p_rhd_pmr has 4 points; 2 and 3 are the central ones
2809	DO_2D( 0, 0, 0, 0)
2810	p_fld_mid(ji,jj,jk) = z_9_16 * ( p_fld_pmr(ji,jj,jk,2) + p_fld_pmr(ji,jj,jk,3) )
2811	END_2D
2812
2813	IF ( k_uv == 1 ) THEN
2814	DO_2D( 0, 0, 0, 0)
2815	IF ( umask(ji-1, jj, jk) > 0.5 ) THEN
2816	z_cor = p_fld_pmr(ji,jj,jk,1)
2817	ELSE
2818	z_cor = p_fld_pmr(ji,jj,jk,2)
2819	END IF
2820	IF ( umask(ji+1, jj, jk) > 0.5 ) THEN
2821	z_cor = z_cor + p_fld_pmr(ji,jj,jk,4)
2822	ELSE
2823	z_cor = z_cor + p_fld_pmr(ji,jj,jk,3)
2824	END IF
2825	p_fld_mid(ji,jj,jk) = p_fld_mid(ji,jj,jk) - z_1_16 * z_cor
2826	END_2D
2827	ELSE ! k_uv == 2
2828	DO_2D( 0, 0, 0, 0)
2829	IF ( vmask(ji, jj-1, jk) > 0.5 ) THEN
2830	z_cor = p_fld_pmr(ji,jj,jk,1)
2831	ELSE
2832	z_cor = p_fld_pmr(ji,jj,jk,2)
2833	END IF
2834	IF ( vmask(ji, jj+1, jk) > 0.5 ) THEN
2835	z_cor = z_cor + p_fld_pmr(ji,jj,jk,4)
2836	ELSE
2837	z_cor = z_cor + p_fld_pmr(ji,jj,jk,3)
2838	END IF
2839	p_fld_mid(ji,jj,jk) = p_fld_mid(ji,jj,jk) - z_1_16 * z_cor
2840	END_2D
2841	END IF ! k_uv
2842
2843	END DO ! jk
2844
2845	END SUBROUTINE calc_mid_cub
2846
2847	!----------------------------------------------------------------------------
2848
2849	SUBROUTINE calc_dz_dij_ccs( k_uv, p_z_pmr, p_dz_ij )
2850
2851	!!---------------------------------------------------------------------
2852	!! * ROUTINE calc_dz_dij_ccs *
2853	!!
2854	!! ** Method : Use constrained cubic spline to determine horizontal derivatives of z at 3 central points
2855	!! The routine is limited to z derivatives because the output is valid on w-levels
2856	!! ** Action : - set p_dz_ij
2857	!!----------------------------------------------------------------------
2858
2859	INTEGER , INTENT(in) :: k_uv ! 1 for u-vel; 2 for v_vel
2860	REAL(wp), DIMENSION(A2D(nn_hls),jpk,4), INTENT(in) :: p_z_pmr ! z field in pmr form
2861	REAL(wp), DIMENSION(A2D(nn_hls),jpk,3), INTENT(out) :: p_dz_ij ! constrained cubic horizontal derivatives at -1/2, 0 and 1/2
2862
2863	REAL(wp), DIMENSION(A2D(nn_hls)) :: zdfld_21, zdfld_32, zdfld_43 ! primitive horizontal differences
2864
2865	INTEGER ji, jj, jk ! standard loop indices
2866	INTEGER jk_msk ! jk level to use for umask and vmask (we need z on the upper and lower faces)
2867	INTEGER j_w_levs ! indicator for data on w-levels
2868
2869	j_w_levs = 2
2870
2871	DO jk = 1, jpk
2872
2873	CALL calc_dfld_pmr_ij( k_uv, j_w_levs, jk, aco_bc_z_hor, bco_bc_z_hor, p_z_pmr, p_dz_ij )
2874
2875	! copy the 2nd horizontal element to the 3rd element of the output array (for an isolated velocity cell p_dz_ij is the same for all 3 elements)
2876	DO_2D( 0, 0, 0, 0)
2877	p_dz_ij(ji,jj,jk,3) = p_dz_ij(ji,jj,jk,2)
2878	END_2D
2879
2880	! set the central element if it is not an isolated velocity cell; this evaluates dz_dij at zeta = 0 ; that is f^{(1)}; and that is given by SMcW(5.8)
2881	IF ( jk == 1 ) THEN
2882	jk_msk = 1
2883	ELSE
2884	jk_msk = jk - 1
2885	END IF
2886
2887	IF ( k_uv == 1 ) THEN
2888	DO_2D( 0, 0, 0, 0)
2889	IF ( umask(ji-1, jj, jk_msk) > 0.5 .OR. umask(ji+1, jj, jk_msk) > 0.5 ) THEN
2890	p_dz_ij(ji,jj,jk,2) = 1.5_wp( p_z_pmr(ji,jj,jk,3) - p_z_pmr(ji,jj,jk,2) ) - 0.25_wp ( p_dz_ij(ji,jj,jk,3) + p_dz_ij(ji,jj,jk,1) )
2891	END IF
2892	END_2D
2893	ELSE ! k_uv == 2
2894	DO_2D( 0, 0, 0, 0)
2895	IF ( vmask(ji, jj-1, jk_msk) > 0.5 .OR. vmask(ji, jj+1, jk_msk) > 0.5 ) THEN
2896	p_dz_ij(ji,jj,jk,2) = 1.5_wp( p_z_pmr(ji,jj,jk,3) - p_z_pmr(ji,jj,jk,2) ) - 0.25_wp ( p_dz_ij(ji,jj,jk,3) + p_dz_ij(ji,jj,jk,1) )
2897	END IF
2898	END_2D
2899	END IF ! k_uv
2900
2901	END DO ! jk
2902
2903	END SUBROUTINE calc_dz_dij_ccs
2904
2905	!----------------------------------------------------------------------------
2906
2907	SUBROUTINE calc_dz_dij_cub( k_uv, p_z_pmr, p_dz_ij )
2908
2909	!!---------------------------------------------------------------------
2910	!! * ROUTINE calc_dz_dij_cub *
2911	!!
2912	!! ** Method : Use simple cubic polynomial to determine horizontal derivatives at 3 central points
2913	!! based on equations (5.5) and (5.6) of SMcW 2003
2914	!! The routine is limited to z derivatives because the output is valid on w-levels
2915	!! ** Action : - set p_dz_ij
2916	!!----------------------------------------------------------------------
2917
2918	INTEGER , INTENT(in) :: k_uv ! 1 for u-vel; 2 for v_vel
2919	REAL(wp), DIMENSION(A2D(nn_hls),jpk,4), INTENT(in) :: p_z_pmr ! z field in pmr form
2920	REAL(wp), DIMENSION(A2D(nn_hls),jpk,3), INTENT(out) :: p_dz_ij ! horizontal derivatives at -1/2, 0 and 1/2
2921
2922	REAL(wp), DIMENSION(A2D(nn_hls)) :: z_z_a, z_z_d ! z fields with boundary conditions applied
2923	REAL(wp) z_1_24 ! 1/24
2924	REAL(wp) z_a, z_b, z_c, z_d ! values of input field at the four points used (-3/2, -1/2, 1/2, 3/2)
2925	REAL(wp) z_c_m_b, z_d_m_a ! differences c - b and d - a
2926	REAL(wp) z_co_1, z_co_2, z_co_3 ! coefficients of polynomial (field = z_co_0 + z_co_1 zeta + .. )
2927	INTEGER ji, jj, jk ! standard loop indices
2928	INTEGER jk_msk ! jk level to use for umask and vmask (we need z on the upper and lower faces)
2929
2930	z_1_24 = 1.0_wp / 24.0_wp
2931
2932	DO jk = 1, jpk
2933
2934	!------------------------------------------------------
2935	! 1. Use simple von Neumann conditions at the boundaries
2936	!------------------------------------------------------
2937	IF ( jk == 1 ) THEN
2938	jk_msk = 1
2939	ELSE
2940	jk_msk = jk - 1
2941	END IF
2942
2943	IF ( k_uv == 1 ) THEN
2944	DO_2D( 0, 0, 0, 0)
2945	IF ( umask(ji-1, jj, jk_msk) > 0.5 ) THEN
2946	z_z_a(ji,jj) = p_z_pmr(ji,jj,jk,1)
2947	ELSE
2948	z_z_a(ji,jj) = p_z_pmr(ji,jj,jk,2)
2949	END IF
2950	IF ( umask(ji+1, jj, jk_msk) > 0.5 ) THEN
2951	z_z_d(ji,jj) = p_z_pmr(ji,jj,jk,4)
2952	ELSE
2953	z_z_d(ji,jj) = p_z_pmr(ji,jj,jk,3)
2954	END IF
2955	END_2D
2956	ELSE ! k_uv == 2
2957	DO_2D( 0, 0, 0, 0)
2958	IF ( vmask(ji, jj-1, jk_msk) > 0.5 ) THEN
2959	z_z_a(ji,jj) = p_z_pmr(ji,jj,jk,1)
2960	ELSE
2961	z_z_a(ji,jj) = p_z_pmr(ji,jj,jk,2)
2962	END IF
2963	IF ( vmask(ji, jj+1, jk_msk) > 0.5 ) THEN
2964	z_z_d(ji,jj) = p_z_pmr(ji,jj,jk,4)
2965	ELSE
2966	z_z_d(ji,jj) = p_z_pmr(ji,jj,jk,3)
2967	END IF
2968	END_2D
2969	END IF ! k_uv
2970
2971	!------------------------------------------------------
2972	! 2. calculate the coefficients for the polynomial fit (c_1, c_2 and c_3; we don't need c_0)
2973	!------------------------------------------------------
2974
2975	DO_2D( 0, 0, 0, 0)
2976	z_a = z_z_a(ji,jj)
2977	z_b = p_z_pmr(ji,jj,jk,2)
2978	z_c = p_z_pmr(ji,jj,jk,3)
2979	z_d = z_z_d(ji,jj)
2980
2981	z_c_m_b = z_c - z_b ; z_d_m_a = z_d - z_a
2982
2983	! eqn (5.6) of SMcW
2984	z_co_1 = 1.125_wp * z_c_m_b - z_1_24 * z_d_m_a
2985	z_co_2 = -0.5_wp * (z_c+z_b) + 0.5_wp * (z_d+z_a)
2986	z_co_3 = -3.0_wp * z_c_m_b + z_d_m_a
2987
2988	! eqn (5.5) of SMcW (first derivative of it)
2989	p_dz_ij(ji,jj,jk,1) = z_co_1 - 0.5_wpz_co_2 + 0.125_wpz_co_3 ! zeta = -0.5
2990	p_dz_ij(ji,jj,jk,2) = z_co_1 ! zeta = 0.0
2991	p_dz_ij(ji,jj,jk,3) = z_co_1 + 0.5_wpz_co_2 + 0.125_wpz_co_3 ! zeta = 0.5
2992
2993	END_2D
2994
2995	END DO ! jk
2996
2997	END SUBROUTINE calc_dz_dij_cub
2998
2999	!----------------------------------------------------------------------------
3000
3001	SUBROUTINE calc_dfld_pmr_ij( k_uv, k_lev_type, kk, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_dfld_ij )
3002
3003	!!---------------------------------------------------------------------
3004	!! * ROUTINE calc_dfld_pmr_ij *
3005	!!
3006	!! ** Method : Calculate constrained spline horizontal derivatives
3007	!!
3008	!! ** Action : - set p_dfld_ij
3009	!!----------------------------------------------------------------------
3010
3011	INTEGER , INTENT(in) :: k_uv ! 1 for u-vel; 2 for v_vel
3012	INTEGER , INTENT(in) :: k_lev_type ! 1 for t-level data; 2 for w-level data; used with check of land/sea mask
3013	INTEGER , INTENT(in) :: kk ! vertical level
3014	REAL(wp) , INTENT(in) :: p_aco_bc_fld_hor ! a coeff for horizontal bc (von Neumann or linear extrapolation)
3015	REAL(wp) , INTENT(in) :: p_bco_bc_fld_hor ! b coeff for horizontal bc (von Neumann or linear extrapolation)
3016	REAL(wp), DIMENSION(A2D(nn_hls),jpk,4), INTENT(in) :: p_fld_pmr ! field in pmr form
3017	REAL(wp), DIMENSION(A2D(nn_hls),jpk,3), INTENT(out) :: p_dfld_ij ! constrained spline horizontal derivatives of the field; only 1 and 2 are set!
3018
3019	REAL(wp), DIMENSION(A2D(nn_hls)) :: zdfld_21, zdfld_32, zdfld_43 ! primitive horizontal differences
3020
3021	INTEGER ji, jj ! standard loop indices
3022	INTEGER jio, jjo ! offsets
3023	INTEGER iktb ! index of the bottom of ref profile
3024
3025	REAL(wp) z1_cff, z_cff_31, z_cff_42 ! temporary sum and products
3026	REAL(wp) zep
3027
3028	INTEGER :: j_t_levs, jk_msk ! indicators that field passed in is valid on t-levels or w-levels; jk to use with masks
3029
3030	j_t_levs = 1
3031
3032	!----------------------------------------------------------------------------------------
3033	! 1. compute and store elementary horizontal differences zfor z_rhd_pmr arrays
3034	!----------------------------------------------------------------------------------------
3035
3036	IF ( k_uv == 1) THEN
3037	jio = 1
3038	jjo = 0
3039	ELSE
3040	jio = 0
3041	jjo = 1
3042	END IF
3043
3044	IF ( k_lev_type == j_t_levs ) THEN
3045	jk_msk = kk
3046	ELSE
3047	IF ( kk == 1 ) THEN
3048	jk_msk = 1
3049	ELSE
3050	jk_msk = kk - 1
3051	END IF
3052	END IF
3053
3054	DO_2D( 0, 0, 0, 0 )
3055	zdfld_21(ji,jj) = p_fld_pmr(ji,jj,kk,2) - p_fld_pmr(ji,jj,kk,1)
3056	zdfld_32(ji,jj) = p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2)
3057	zdfld_43(ji,jj) = p_fld_pmr(ji,jj,kk,4) - p_fld_pmr(ji,jj,kk,3)
3058	END_2D
3059
3060	zep = 1.e-15
3061	DO_2D( 0, 0, 0, 0 )
3062	z_cff_31 = MAX( 2._wp * zdfld_21(ji,jj) * zdfld_32(ji,jj), 0._wp )
3063	z1_cff = zdfld_21(ji,jj) + zdfld_32(ji,jj)
3064	p_dfld_ij(ji,jj,kk,1) = z_cff_31 / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
3065
3066	z_cff_42 = MAX( 2._wp * zdfld_32(ji,jj) * zdfld_43(ji,jj), 0._wp )
3067	z1_cff = zdfld_32(ji,jj) + zdfld_43(ji,jj)
3068	p_dfld_ij(ji,jj,kk,2) = z_cff_42 / SIGN( MAX( ABS(z1_cff), zep ), z1_cff )
3069	END_2D
3070
3071	!----------------------------------------------------------------------------------
3072	! 2. apply boundary conditions at side boundaries using 5.36-5.37
3073	!----------------------------------------------------------------------------------
3074
3075	IF ( k_uv == 1 ) THEN
3076	DO_2D( 0, 0, 0, 0)
3077	! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj)
3078	IF ( umask(ji,jj,jk_msk) > 0.5_wp .AND. umask(ji-1,jj,jk_msk) < 0.5_wp .AND. umask(ji+1,jj,jk_msk) > 0.5_wp) THEN
3079	p_dfld_ij(ji,jj,kk,1) = p_aco_bc_fld_hor * ( p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2) ) - p_bco_bc_fld_hor * p_dfld_ij(ji,jj,kk,2)
3080	END IF
3081	! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj)
3082	IF ( umask(ji,jj,jk_msk) < 0.5_wp .AND. umask(ji-1,jj,jk_msk) > 0.5_wp .AND. umask(ji-2,jj,jk_msk) > 0.5_wp) THEN
3083	p_dfld_ij(ji,jj,kk,2) = p_aco_bc_fld_hor * ( p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2) ) - p_bco_bc_fld_hor * p_dfld_ij(ji,jj,kk,1)
3084	END IF
3085	END_2D
3086
3087	DO_2D( 0, 0, 0, 0)
3088	! For an isolated velocity point use the central difference (this is important)
3089	IF ( umask(ji-1, jj, jk_msk) < 0.5 .AND. umask(ji+1, jj, jk_msk) < 0.5 ) THEN
3090	p_dfld_ij(ji,jj,kk,1) = p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2)
3091	p_dfld_ij(ji,jj,kk,2) = p_dfld_ij(ji,jj,kk,1)
3092	END IF
3093	END_2D
3094
3095	ELSE ! k_uv == 2
3096
3097	DO_2D( 0, 0, 0, 0)
3098	! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi)
3099	IF ( vmask(ji,jj,jk_msk) > 0.5_wp .AND. vmask(ji,jj-1,jk_msk) < 0.5_wp .AND. vmask(ji,jj+1,jk_msk) > 0.5_wp) THEN
3100	p_dfld_ij(ji,jj,kk,1) = p_aco_bc_fld_hor * (p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2) ) - p_bco_bc_fld_hor * p_dfld_ij(ji,jj,kk,2)
3101	END IF
3102	! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi)
3103	IF ( vmask(ji,jj,jk_msk) < 0.5_wp .AND. vmask(ji,jj-1,jk_msk) > 0.5_wp .AND. vmask(ji,jj-2,jk_msk) > 0.5_wp) THEN
3104	p_dfld_ij(ji,jj,kk,2) = p_aco_bc_fld_hor * (p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2) ) - p_bco_bc_fld_hor * p_dfld_ij(ji,jj,kk,1)
3105	END IF
3106	END_2D
3107
3108	DO_2D( 0, 0, 0, 0)
3109	! For an isolated velocity point use the central difference (this is important)
3110	IF ( vmask(ji, jj-1, jk_msk) < 0.5 .AND. vmask(ji, jj+1, jk_msk) < 0.5 ) THEN
3111	p_dfld_ij(ji,jj,kk,1) = p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2)
3112	p_dfld_ij(ji,jj,kk,2) = p_dfld_ij(ji,jj,kk,1)
3113	END IF
3114	END_2D
3115
3116	END IF ! k_uv
3117
3118	END SUBROUTINE calc_dfld_pmr_ij
3119
3120	!------------------------------------------------------------------------------------------
3121
3122	SUBROUTINE vrt_int_quad(k_mbkt, p_rhd, p_e3t, p_p, p_F)
3123
3124	!!---------------------------------------------------------------------
3125	!! * ROUTINE calc_rhd_k *
3126	!!
3127	!! ** Method : Use quadratic reconstruction of p_rhd (density profile) to calculate pressure profile and
3128	!! forces on vertical faces between w levels
3129	!!
3130	!! ** Action : - set p_p and p_F
3131	!!----------------------------------------------------------------------
3132
3133	INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: k_mbkt ! bottom level
3134	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_rhd ! densities on tracer levels
3135	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_e3t ! layer thickness between w levels
3136	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(out) :: p_p ! pressures on w levels
3137	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(out) :: p_F ! force on the vertical face between w levels
3138
3139	INTEGER :: ji, jj, jk ! dummy loop indices
3140	REAL(wp) :: z_2_3, z_4_3, z_8_3 ! 2/3 and 4/3 and 8/3
3141	REAL(wp) :: z_22_3, z_28_3 ! 22/3 and 28/3
3142
3143	REAL(wp) :: zr_a, zr_b, zr_c ! rhd evaluated at i, i+1/2 and i+1
3144	REAL(wp) :: za_0, za_1, za_2 ! polynomial coefficients
3145	REAL(wp) :: ze3t !
3146
3147	z_2_3 = 2.0_wp/3.0_wp
3148	z_4_3 = 4.0_wp/3.0_wp
3149	z_8_3 = 8.0_wp/3.0_wp
3150	z_22_3 = 22.0_wp/3.0_wp
3151	z_28_3 = 28.0_wp/3.0_wp
3152
3153	! Integrate densities in the vertical using quadratic fits to the densities
3154	! zr_a is r_{-1} zr_b is r_1 ; zr_c is r_3
3155
3156	! one-sided quadratic at the upper boundary. pressure is zero at the upper surface
3157	DO_2D(0,0,0,0)
3158	zr_a = p_rhd(ji,jj,1)
3159	zr_b = p_rhd(ji,jj,2)
3160	zr_c = p_rhd(ji,jj,3)
3161	ze3t = p_e3t(ji,jj,1)
3162	za_0 = 0.125_wp * ( - zr_c + 6._wpzr_b + 3._wpzr_a )
3163	za_1 = 0.5_wp * ( zr_b - zr_a )
3164	za_2 = 0.125_wp * ( zr_c - 2._wp*zr_b + zr_a )
3165	p_p(ji,jj,1) = 0.0_wp
3166	p_p(ji,jj,2) = gravze3t ( za_0 - za_1 + z_4_3*za_2 )
3167	p_F(ji,jj,1) = 0.5_wpgravze3tze3t( za_0 - z_4_3za_1 + 2._wpza_2 )
3168	END_2D
3169
3170	! centred quadratic in the interior
3171	DO_3D(0,0,0,0,2,jpk-1)
3172	zr_a = p_rhd(ji,jj,jk-1)
3173	zr_b = p_rhd(ji,jj,jk)
3174	zr_c = p_rhd(ji,jj,jk+1)
3175	ze3t = p_e3t(ji,jj,jk)
3176	za_0 = 0.125_wp * ( - zr_c + 6._wpzr_b + 3._wpzr_a )
3177	za_1 = 0.5_wp * ( zr_b - zr_a )
3178	za_2 = 0.125_wp * ( zr_c - 2._wp*zr_b + zr_a )
3179	p_p(ji,jj,jk+1) = p_p(ji,jj,jk) + gravze3t ( za_0 + za_1 + z_4_3*za_2 )
3180	p_F(ji,jj,jk ) = ze3tp_p(ji,jj,jk) + 0.5gravze3tze3t( za_0 + z_2_3za_1 + z_2_3*za_2 )
3181	END_3D
3182
3183	! one-sided quadratic at the lower boundary
3184	DO_2D(0,0,0,0)
3185	jk = k_mbkt(ji,jj)
3186	zr_a = p_rhd(ji,jj,jk-2)
3187	zr_b = p_rhd(ji,jj,jk-1)
3188	zr_c = p_rhd(ji,jj,jk)
3189	ze3t = p_e3t(ji,jj,jk)
3190	za_0 = 0.125_wp * ( - zr_c + 6._wpzr_b + 3._wpzr_a )
3191	za_1 = 0.5_wp * ( zr_b - zr_a )
3192	za_2 = 0.125_wp * ( zr_c - 2._wp*zr_b + zr_a )
3193	p_p(ji,jj,jk+1) = p_p(ji,jj,jk) + gravze3t ( za_0 + 3.za_1 + z_28_3za_2 )
3194	p_F(ji,jj,jk ) = ze3tp_p(ji,jj,jk) + 0.5gravze3tze3t( za_0 + z_8_3za_1 + z_22_3*za_2 )
3195	END_2D
3196
3197	RETURN
3198	END SUBROUTINE vrt_int_quad
3199
3200	!------------------------------------------------------------------------------------------
3201
3202	SUBROUTINE dbg_2di( cc_array, ki_2d )
3203	CHARACTER(), INTENT(IN) :: cc_array
3204	INTEGER, DIMENSION(A2D(nn_hls)), INTENT(IN) :: ki_2d
3205
3206	INTEGER :: ji_prt_low, ji_prt_upp, jj_prt_low, jj_prt_upp
3207	INTEGER :: ji_prt, jj_prt
3208
3209	IF( lwp ) THEN
3210	ji_prt_low = MAX( ntsi-nn_hls, ki_dbg_min )
3211	ji_prt_upp = MIN( ntei+nn_hls, ki_dbg_max )
3212	jj_prt_low = MAX( ntsj-nn_hls, kj_dbg_min )
3213	jj_prt_upp = MIN( ntej+nn_hls, kj_dbg_max )
3214
3215	! print out a 2D horizontal slice
3216
3217	IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ij ) THEN
3218	IF ( jj_prt_low <= jj_prt_upp ) THEN
3219	WRITE(numout,*)
3220	WRITE(numout,*) cc_array
3221	WRITE(numout,*) ' ji_prt_low, ji_prt_upp = ', ji_prt_low, ji_prt_upp
3222	WRITE(numout,*) ' row/col ', ( ji_prt, ji_prt = ji_prt_low, ji_prt_upp )
3223	DO jj_prt = jj_prt_low, jj_prt_upp
3224	WRITE(numout,*) jj_prt, ( ki_2d(ji_prt, jj_prt), ji_prt = ji_prt_low, ji_prt_upp )
3225	END DO
3226	END IF
3227	END IF
3228
3229	END IF
3230
3231	RETURN
3232	END SUBROUTINE dbg_2di
3233
3234	!----------------------------------------------------------------------------
3235
3236	SUBROUTINE dbg_2di_k( cc_array, ki_2d, kk )
3237	CHARACTER(), INTENT(IN) :: cc_array
3238	INTEGER, DIMENSION(A2D(nn_hls)), INTENT(IN) :: ki_2d
3239	INTEGER, INTENT(IN) :: kk
3240
3241	INTEGER :: ji_prt_low, ji_prt_upp, jj_prt_low, jj_prt_upp
3242	INTEGER :: ji_prt, jj_prt, jk_prt
3243
3244	IF( lwp ) THEN
3245	ji_prt_low = MAX( ntsi, ki_dbg_min )
3246	ji_prt_upp = MIN( ntei, ki_dbg_max )
3247	jj_prt_low = MAX( ntsj, kj_dbg_min )
3248	jj_prt_upp = MIN( ntej, kj_dbg_max )
3249
3250	! print out a 2D (ji, jk) slice
3251
3252	IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ik ) THEN
3253	IF ( ntsj <= kj_dbg_cen .AND. kj_dbg_cen <= ntej ) THEN
3254	IF ( kk == kk_dbg_min ) THEN
3255	WRITE(numout,*)
3256	WRITE(numout,*) cc_array, ' kj = ', kj_dbg_cen
3257	WRITE(numout,*) 'ji_prt_low, ji_prt_upp = ', ji_prt_low, ji_prt_upp
3258	WRITE(numout,*) 'kk_dbg_min, kk_dbg_max = ', kk_dbg_min, kk_dbg_max
3259	END IF
3260	IF ( kk_dbg_min <= kk .AND. kk <= kk_dbg_max ) &
3261	& WRITE(numout,*) cc_array, kk, ( ki_2d(ji_prt, kj_dbg_cen), ji_prt = ji_prt_low, ji_prt_upp )
3262	END IF
3263	END IF
3264
3265	! print out a 2D (jj, jk) slice
3266
3267	IF ( jj_prt_low <= jj_prt_upp .AND. ln_dbg_jk ) THEN
3268	IF ( ntsi <= ki_dbg_cen .AND. ki_dbg_cen <= ntei ) THEN
3269	IF ( kk == kk_dbg_min ) THEN
3270	WRITE(numout,*)
3271	WRITE(numout,*) cc_array, ' ki = ', ki_dbg_cen
3272	WRITE(numout,*) 'jj_prt_low, jj_prt_upp = ', jj_prt_low, jj_prt_upp
3273	WRITE(numout,*) 'kk_dbg_min, kk_dbg_max = ', kk_dbg_min, kk_dbg_max
3274	END IF
3275	IF ( kk_dbg_min <= kk .AND. kk <= kk_dbg_max ) &
3276	& WRITE(numout,*) cc_array, kk, ( ki_2d(ki_dbg_cen, jj_prt), jj_prt = jj_prt_low, jj_prt_upp )
3277	END IF
3278	END IF
3279
3280	END IF
3281
3282	RETURN
3283	END SUBROUTINE dbg_2di_k
3284
3285	!----------------------------------------------------------------------------
3286
3287	SUBROUTINE dbg_2dr( cc_array, pr_2d )
3288	CHARACTER(), INTENT(IN) :: cc_array
3289	REAL(wp), DIMENSION(A2D(nn_hls)), INTENT(IN) :: pr_2d
3290
3291	INTEGER :: ji_prt_low, ji_prt_upp, jj_prt_low, jj_prt_upp
3292	INTEGER :: ji_prt, jj_prt
3293
3294	IF(lwp) THEN
3295	ji_prt_low = MAX( ntsi-nn_hls, ki_dbg_min )
3296	ji_prt_upp = MIN( ntei+nn_hls, ki_dbg_max )
3297	jj_prt_low = MAX( ntsj-nn_hls, kj_dbg_min )
3298	jj_prt_upp = MIN( ntej+nn_hls, kj_dbg_max )
3299
3300	! print out a 2D horizontal slice
3301
3302	IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ij ) THEN
3303	IF ( jj_prt_low <= jj_prt_upp ) THEN
3304	WRITE(numout,*)
3305	WRITE(numout,*) cc_array
3306	WRITE(numout,*) ' row/col ', ( ji_prt, ji_prt = ji_prt_low, ji_prt_upp )
3307	DO jj_prt = jj_prt_low, jj_prt_upp
3308	WRITE(numout,*) jj_prt, ( pr_2d(ji_prt, jj_prt), ji_prt = ji_prt_low, ji_prt_upp )
3309	END DO
3310	END IF
3311	END IF
3312
3313	END IF
3314
3315	RETURN
3316	END SUBROUTINE dbg_2dr
3317
3318	!----------------------------------------------------------------------------
3319
3320	SUBROUTINE dbg_3di( cc_array, ki_3d )
3321
3322	CHARACTER(), INTENT(IN) :: cc_array
3323	INTEGER, DIMENSION(A2D(nn_hls),jpk), INTENT(IN) :: ki_3d
3324
3325	INTEGER :: ji_prt_low, ji_prt_upp, jj_prt_low, jj_prt_upp
3326	INTEGER :: ji_prt, jj_prt, jk_prt
3327
3328	! output values
3329
3330
3331	IF(lwp) THEN
3332	ji_prt_low = MAX( ntsi-nn_hls, ki_dbg_min )
3333	ji_prt_upp = MIN( ntei+nn_hls, ki_dbg_max )
3334	jj_prt_low = MAX( ntsj-nn_hls, kj_dbg_min )
3335	jj_prt_upp = MIN( ntej+nn_hls, kj_dbg_max )
3336
3337	! print out a 2D horizontal slice
3338
3339	IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ij ) THEN
3340	IF ( jj_prt_low <= jj_prt_upp ) THEN
3341	WRITE(numout,*)
3342	WRITE(numout,*) cc_array, ' level = ', kk_dbg_cen
3343	WRITE(numout,*) ' row/col ', ( ji_prt, ji_prt = ji_prt_low, ji_prt_upp )
3344	DO jj_prt = jj_prt_low, jj_prt_upp
3345	WRITE(numout,*) jj_prt, ( ki_3d(ji_prt, jj_prt, kk_dbg_cen), ji_prt = ji_prt_low, ji_prt_upp )
3346	END DO
3347	END IF
3348	END IF
3349
3350	! print out a 2D (ji, jk) slice
3351
3352	IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ik ) THEN
3353	IF ( ntsj <= kj_dbg_cen .AND. kj_dbg_cen <= ntej ) THEN
3354	WRITE(numout,*)
3355	WRITE(numout,*) cc_array, ' kj = ', kj_dbg_cen
3356	WRITE(numout,*) ' row/lev ', ( ji_prt, ji_prt = ji_prt_low, ji_prt_upp )
3357	DO jk_prt = kk_dbg_min, kk_dbg_max
3358	WRITE(numout,*) jk_prt, ( ki_3d(ji_prt, kj_dbg_cen, jk_prt), ji_prt = ji_prt_low, ji_prt_upp )
3359	END DO
3360	END IF
3361	END IF
3362
3363	! print out a 2D (jj, jk) slice
3364
3365	IF ( jj_prt_low <= jj_prt_upp .AND. ln_dbg_jk ) THEN
3366	IF ( ntsi <= ki_dbg_cen .AND. ki_dbg_cen <= ntei ) THEN
3367	WRITE(numout,*)
3368	WRITE(numout,*) cc_array, ' ki = ', ki_dbg_cen
3369	WRITE(numout,*) ' col/lev ', ( jj_prt, jj_prt = jj_prt_low, jj_prt_upp )
3370	DO jk_prt = kk_dbg_min, kk_dbg_max
3371	WRITE(numout,*) jk_prt, ( ki_3d(ki_dbg_cen, jj_prt, jk_prt), jj_prt = jj_prt_low, jj_prt_upp )
3372	END DO
3373	END IF
3374	END IF
3375
3376	END IF
3377
3378	RETURN
3379	END SUBROUTINE dbg_3di
3380
3381	!----------------------------------------------------------------------------
3382
3383	SUBROUTINE dbg_3dr( cc_array, pr_3d )
3384
3385	CHARACTER(), INTENT(IN) :: cc_array
3386	REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(IN) :: pr_3d
3387
3388	INTEGER :: ji_prt_low, ji_prt_upp, jj_prt_low, jj_prt_upp
3389	INTEGER :: ji_prt, jj_prt, jk_prt
3390
3391	! output values
3392
3393
3394	IF(lwp) THEN
3395	ji_prt_low = MAX( ntsi-nn_hls, ki_dbg_min )
3396	ji_prt_upp = MIN( ntei+nn_hls, ki_dbg_max )
3397	jj_prt_low = MAX( ntsj-nn_hls, kj_dbg_min )
3398	jj_prt_upp = MIN( ntej+nn_hls, kj_dbg_max )
3399
3400	! print out a 2D horizontal slice
3401
3402	IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ij ) THEN
3403	IF ( jj_prt_low <= jj_prt_upp ) THEN
3404	WRITE(numout,*)
3405	WRITE(numout,*) cc_array, ' level = ', kk_dbg_cen
3406	WRITE(numout,*) ' row/col ', ( ji_prt, ji_prt = ji_prt_low, ji_prt_upp )
3407	DO jj_prt = jj_prt_low, jj_prt_upp
3408	WRITE(numout,*) jj_prt, ( pr_3d(ji_prt, jj_prt, kk_dbg_cen), ji_prt = ji_prt_low, ji_prt_upp )
3409	END DO
3410	END IF
3411	END IF
3412
3413	! print out a 2D (ji, jk) slice
3414
3415	IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ik ) THEN
3416	IF ( ntsj <= kj_dbg_cen .AND. kj_dbg_cen <= ntej ) THEN
3417	WRITE(numout,*)
3418	WRITE(numout,*) cc_array, ' kj = ', kj_dbg_cen
3419	WRITE(numout,*) ' row/lev ', ( ji_prt, ji_prt = ji_prt_low, ji_prt_upp )
3420	DO jk_prt = kk_dbg_min, kk_dbg_max
3421	WRITE(numout,*) jk_prt, ( pr_3d(ji_prt, kj_dbg_cen, jk_prt), ji_prt = ji_prt_low, ji_prt_upp )
3422	END DO
3423	END IF
3424	END IF
3425
3426	! print out a 2D (jj, jk) slice
3427
3428	IF ( jj_prt_low <= jj_prt_upp .AND. ln_dbg_jk ) THEN
3429	IF ( ntsi <= ki_dbg_cen .AND. ki_dbg_cen <= ntei ) THEN
3430	WRITE(numout,*)
3431	WRITE(numout,*) cc_array, ' ki = ', ki_dbg_cen
3432	WRITE(numout,*) ' col/lev ', ( jj_prt, jj_prt = jj_prt_low, jj_prt_upp )
3433	DO jk_prt = kk_dbg_min, kk_dbg_max
3434	WRITE(numout,*) jk_prt, ( pr_3d(ki_dbg_cen, jj_prt, jk_prt), jj_prt = jj_prt_low, jj_prt_upp )
3435	END DO
3436	END IF
3437	END IF
3438
3439	END IF
3440
3441	RETURN
3442	END SUBROUTINE dbg_3dr
3443
3444
3445	!!======================================================================
3446	END MODULE dynhpg

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

Download in other formats: