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

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

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

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

Download in other formats: