Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

dynhpg.F90 @ 15035

Last change on this file since 15035 was 15035, checked in by ayoung, 3 years ago
Backporting djc into r4.0-HEAD. #2480.
Property svn:keywords set to `Id`
File size: 79.2 KB

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

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

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

Context Navigation

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

Download in other formats: