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 @ 15728

Last change on this file since 15728 was 15728, checked in by dbruciaferri, 2 years ago
updating to V3 of DJR and FFR code - dynhpg_djr_ffr_smooth_v3.F90
Property svn:keywords set to `Id`
File size: 196.3 KB

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

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

Download in other formats: