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 branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2015/dev_r5803_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 5870

Last change on this file since 5870 was 5870, checked in by acc, 8 years ago
Branch 2015/dev_r5803_NOC_WAD. Merge in trunk changes from 5803 to 5869 in preparation for merge. Also tidied and reorganised some wetting and drying code. Renamed wadlmt.F90 to wetdry.F90. Wetting drying code changes restricted to domzgr.F90, domvvl.F90 nemogcm.F90 sshwzv.F90, dynspg_ts.F90, wetdry.F90 and dynhpg.F90. Code passes full SETTE tests with ln_wd=.false.. Still awaiting test case for checking with ln_wd=.false.
Property svn:keywords set to `Id`
File size: 77.8 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	!! 3.6? ! 2015-11 (H. Liu) add Wetting/Drying pressure filter
20	!!----------------------------------------------------------------------
21
22	!!----------------------------------------------------------------------
23	!! dyn_hpg : update the momentum trend with the now horizontal
24	!! gradient of the hydrostatic pressure
25	!! dyn_hpg_init : initialisation and control of options
26	!! hpg_zco : z-coordinate scheme
27	!! hpg_zps : z-coordinate plus partial steps (interpolation)
28	!! hpg_sco : s-coordinate (standard jacobian formulation)
29	!! hpg_isf : s-coordinate (sco formulation) adapted to ice shelf
30	!! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial)
31	!! hpg_prj : s-coordinate (Pressure Jacobian with Cubic polynomial)
32	!!----------------------------------------------------------------------
33	USE oce ! ocean dynamics and tracers
34	USE sbc_oce ! surface variable (only for the flag with ice shelf)
35	USE dom_oce ! ocean space and time domain
36	USE wet_dry ! wetting and drying
37	USE phycst ! physical constants
38	USE trd_oce ! trends: ocean variables
39	USE trddyn ! trend manager: dynamics
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 wrk_nemo ! Memory Allocation
47	USE timing ! Timing
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	LOGICAL , PUBLIC :: ln_dynhpg_imp !: semi-implicite hpg flag
63
64	INTEGER , PUBLIC :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM)
65
66	!! * Substitutions
67	# include "domzgr_substitute.h90"
68	# include "vectopt_loop_substitute.h90"
69	!!----------------------------------------------------------------------
70	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
71	!! $Id$
72	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
73	!!----------------------------------------------------------------------
74	CONTAINS
75
76	SUBROUTINE dyn_hpg( kt )
77	!!---------------------------------------------------------------------
78	!! * ROUTINE dyn_hpg *
79	!!
80	!! ** Method : Call the hydrostatic pressure gradient routine
81	!! using the scheme defined in the namelist
82	!!
83	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
84	!! - send trends to trd_dyn for futher diagnostics (l_trddyn=T)
85	!!----------------------------------------------------------------------
86	INTEGER, INTENT(in) :: kt ! ocean time-step index
87	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv
88	!!----------------------------------------------------------------------
89	!
90	IF( nn_timing == 1 ) CALL timing_start('dyn_hpg')
91	!
92	IF( l_trddyn ) THEN ! Temporary saving of ua and va trends (l_trddyn)
93	CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv )
94	ztrdu(:,:,:) = ua(:,:,:)
95	ztrdv(:,:,:) = va(:,:,:)
96	ENDIF
97	!
98	SELECT CASE ( nhpg ) ! Hydrostatic pressure gradient computation
99	CASE ( 0 ) ; CALL hpg_zco ( kt ) ! z-coordinate
100	CASE ( 1 ) ; CALL hpg_zps ( kt ) ! z-coordinate plus partial steps (interpolation)
101	CASE ( 2 ) ; CALL hpg_sco ( kt ) ! s-coordinate (standard jacobian formulation)
102	CASE ( 3 ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial)
103	CASE ( 4 ) ; CALL hpg_prj ( kt ) ! s-coordinate (Pressure Jacobian scheme)
104	CASE ( 5 ) ; CALL hpg_isf ( kt ) ! s-coordinate similar to sco modify for ice shelf
105	END SELECT
106	!
107	IF( l_trddyn ) THEN ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
108	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
109	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
110	CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt )
111	CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv )
112	ENDIF
113	!
114	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg - Ua: ', mask1=umask, &
115	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
116	!
117	IF( nn_timing == 1 ) CALL timing_stop('dyn_hpg')
118	!
119	END SUBROUTINE dyn_hpg
120
121
122	SUBROUTINE dyn_hpg_init
123	!!----------------------------------------------------------------------
124	!! * ROUTINE dyn_hpg_init *
125	!!
126	!! ** Purpose : initializations for the hydrostatic pressure gradient
127	!! computation and consistency control
128	!!
129	!! ** Action : Read the namelist namdyn_hpg and check the consistency
130	!! with the type of vertical coordinate used (zco, zps, sco)
131	!!----------------------------------------------------------------------
132	INTEGER :: ioptio = 0 ! temporary integer
133	INTEGER :: ios ! Local integer output status for namelist read
134	!!
135	NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, &
136	& ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, ln_dynhpg_imp
137	!!----------------------------------------------------------------------
138	!
139	REWIND( numnam_ref ) ! Namelist namdyn_hpg in reference namelist : Hydrostatic pressure gradient
140	READ ( numnam_ref, namdyn_hpg, IOSTAT = ios, ERR = 901)
141	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in reference namelist', lwp )
142
143	REWIND( numnam_cfg ) ! Namelist namdyn_hpg in configuration namelist : Hydrostatic pressure gradient
144	READ ( numnam_cfg, namdyn_hpg, IOSTAT = ios, ERR = 902 )
145	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_hpg in configuration namelist', lwp )
146	IF(lwm) WRITE ( numond, namdyn_hpg )
147	!
148	IF(lwp) THEN ! Control print
149	WRITE(numout,*)
150	WRITE(numout,*) 'dyn_hpg_init : hydrostatic pressure gradient initialisation'
151	WRITE(numout,*) '~~~~~~~~~~~~'
152	WRITE(numout,*) ' Namelist namdyn_hpg : choice of hpg scheme'
153	WRITE(numout,*) ' z-coord. - full steps ln_hpg_zco = ', ln_hpg_zco
154	WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps
155	WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco
156	WRITE(numout,*) ' s-coord. (standard jacobian formulation) for isf ln_hpg_isf = ', ln_hpg_isf
157	WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc
158	WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj
159	WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp
160	ENDIF
161	!
162	IF( ln_hpg_djc ) &
163	& CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method &
164	& currently disabled (bugs under investigation). Please select &
165	& either ln_hpg_sco or ln_hpg_prj instead')
166	!
167	IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) &
168	& CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:&
169	& the standard jacobian formulation hpg_sco or &
170	& the pressure jacobian formulation hpg_prj')
171
172	IF( ln_hpg_isf .AND. .NOT. ln_isfcav ) &
173	& CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' )
174	IF( .NOT. ln_hpg_isf .AND. ln_isfcav ) &
175	& CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' )
176	!
177	! ! Set nhpg from ln_hpg_... flags
178	IF( ln_hpg_zco ) nhpg = 0
179	IF( ln_hpg_zps ) nhpg = 1
180	IF( ln_hpg_sco ) nhpg = 2
181	IF( ln_hpg_djc ) nhpg = 3
182	IF( ln_hpg_prj ) nhpg = 4
183	IF( ln_hpg_isf ) nhpg = 5
184	!
185	! ! Consistency check
186	ioptio = 0
187	IF( ln_hpg_zco ) ioptio = ioptio + 1
188	IF( ln_hpg_zps ) ioptio = ioptio + 1
189	IF( ln_hpg_sco ) ioptio = ioptio + 1
190	IF( ln_hpg_djc ) ioptio = ioptio + 1
191	IF( ln_hpg_prj ) ioptio = ioptio + 1
192	IF( ln_hpg_isf ) ioptio = ioptio + 1
193	IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
194	!
195	! initialisation of ice load
196	riceload(:,:)=0.0
197	!
198	END SUBROUTINE dyn_hpg_init
199
200
201	SUBROUTINE hpg_zco( kt )
202	!!---------------------------------------------------------------------
203	!! * ROUTINE hpg_zco *
204	!!
205	!! ** Method : z-coordinate case, levels are horizontal surfaces.
206	!! The now hydrostatic pressure gradient at a given level, jk,
207	!! is computed by taking the vertical integral of the in-situ
208	!! density gradient along the model level from the suface to that
209	!! level: zhpi = grav .....
210	!! zhpj = grav .....
211	!! add it to the general momentum trend (ua,va).
212	!! ua = ua - 1/e1u * zhpi
213	!! va = va - 1/e2v * zhpj
214	!!
215	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
216	!!----------------------------------------------------------------------
217	INTEGER, INTENT(in) :: kt ! ocean time-step index
218	!!
219	INTEGER :: ji, jj, jk ! dummy loop indices
220	REAL(wp) :: zcoef0, zcoef1 ! temporary scalars
221	REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj
222	!!----------------------------------------------------------------------
223	!
224	CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
225	!
226	IF( kt == nit000 ) THEN
227	IF(lwp) WRITE(numout,*)
228	IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend'
229	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case '
230	ENDIF
231
232	zcoef0 = - grav * 0.5_wp ! Local constant initialization
233
234	! Surface value
235	DO jj = 2, jpjm1
236	DO ji = fs_2, fs_jpim1 ! vector opt.
237	zcoef1 = zcoef0 * fse3w(ji,jj,1)
238	! hydrostatic pressure gradient
239	zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
240	zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
241	! add to the general momentum trend
242	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
243	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
244	END DO
245	END DO
246
247	!
248	! interior value (2=<jk=<jpkm1)
249	DO jk = 2, jpkm1
250	DO jj = 2, jpjm1
251	DO ji = fs_2, fs_jpim1 ! vector opt.
252	zcoef1 = zcoef0 * fse3w(ji,jj,jk)
253	! hydrostatic pressure gradient
254	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &
255	& + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) &
256	& - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) / e1u(ji,jj)
257
258	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &
259	& + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) &
260	& - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) / e2v(ji,jj)
261	! add to the general momentum trend
262	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
263	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
264	END DO
265	END DO
266	END DO
267	!
268	CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
269	!
270	END SUBROUTINE hpg_zco
271
272
273	SUBROUTINE hpg_zps( kt )
274	!!---------------------------------------------------------------------
275	!! * ROUTINE hpg_zps *
276	!!
277	!! ** Method : z-coordinate plus partial steps case. blahblah...
278	!!
279	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
280	!!----------------------------------------------------------------------
281	INTEGER, INTENT(in) :: kt ! ocean time-step index
282	!!
283	INTEGER :: ji, jj, jk ! dummy loop indices
284	INTEGER :: iku, ikv ! temporary integers
285	REAL(wp) :: zcoef0, zcoef1, zcoef2, zcoef3 ! temporary scalars
286	REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj
287	!!----------------------------------------------------------------------
288	!
289	CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
290	!
291	IF( kt == nit000 ) THEN
292	IF(lwp) WRITE(numout,*)
293	IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend'
294	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate with partial steps - vector optimization'
295	ENDIF
296
297
298	! Local constant initialization
299	zcoef0 = - grav * 0.5_wp
300
301	! Surface value (also valid in partial step case)
302	DO jj = 2, jpjm1
303	DO ji = fs_2, fs_jpim1 ! vector opt.
304	zcoef1 = zcoef0 * fse3w(ji,jj,1)
305	! hydrostatic pressure gradient
306	zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
307	zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) / 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 * fse3w(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) ) ) / 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) ) ) / 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
336	! partial steps correction at the last level (use gru & grv computed in zpshde.F90)
337	DO jj = 2, jpjm1
338	DO ji = 2, jpim1
339	iku = mbku(ji,jj)
340	ikv = mbkv(ji,jj)
341	zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj ,iku) )
342	zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji ,jj+1,ikv) )
343	IF( iku > 1 ) THEN ! on i-direction (level 2 or more)
344	ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value
345	zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one
346	& + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj)
347	ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend
348	ENDIF
349	IF( ikv > 1 ) THEN ! on j-direction (level 2 or more)
350	va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value
351	zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one
352	& + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj)
353	va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend
354	ENDIF
355	END DO
356	END DO
357	!
358	CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
359	!
360	END SUBROUTINE hpg_zps
361
362	SUBROUTINE hpg_sco( kt )
363	!!---------------------------------------------------------------------
364	!! * ROUTINE hpg_sco *
365	!!
366	!! ** Method : s-coordinate case. Jacobian scheme.
367	!! The now hydrostatic pressure gradient at a given level, jk,
368	!! is computed by taking the vertical integral of the in-situ
369	!! density gradient along the model level from the suface to that
370	!! level. s-coordinates (ln_sco): a corrective term is added
371	!! to the horizontal pressure gradient :
372	!! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ]
373	!! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ]
374	!! add it to the general momentum trend (ua,va).
375	!! ua = ua - 1/e1u * zhpi
376	!! va = va - 1/e2v * zhpj
377	!!
378	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
379	!!----------------------------------------------------------------------
380	INTEGER, INTENT(in) :: kt ! ocean time-step index
381	!!
382	INTEGER :: ji, jj, jk, jii, jjj ! dummy loop indices
383	REAL(wp) :: zcoef0, zuap, zvap, znad, ztmp ! temporary scalars
384	LOGICAL :: ll_tmp1, ll_tmp2, ll_tmp3 ! local logical variables
385	REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj
386	REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter
387	!!----------------------------------------------------------------------
388	!
389	CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj )
390	IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
391	!
392	IF( kt == nit000 ) THEN
393	IF(lwp) WRITE(numout,*)
394	IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
395	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used'
396	ENDIF
397
398	! Local constant initialization
399	zcoef0 = - grav * 0.5_wp
400	! To use density and not density anomaly
401	IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume
402	ELSE ; znad = 0._wp ! Fixed volume
403	ENDIF
404
405	IF(ln_wd) THEN
406	DO jj = 2, jpjm1
407	DO ji = 2, jpim1
408	ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj))
409	ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) > rn_wdmin1 + rn_wdmin2
410	ll_tmp3 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + &
411	& rn_wdmin1 + rn_wdmin2
412
413	IF(ll_tmp1.AND.ll_tmp2) THEN
414	zcpx(ji,jj) = 1.0_wp
415	wduflt(ji,jj) = 1.0_wp
416	ELSE IF(ll_tmp3) THEN
417	! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here
418	zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) / &
419	& (sshn(ji+1,jj) - sshn(ji,jj)))
420	wduflt(ji,jj) = 1.0_wp
421	ELSE
422	zcpx(ji,jj) = 0._wp
423	wduflt(ji,jj) = 0.0_wp
424	END IF
425
426	ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1))
427	ll_tmp2 = MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) > rn_wdmin1 + rn_wdmin2
428	ll_tmp3 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + &
429	& rn_wdmin1 + rn_wdmin2
430
431	IF(ll_tmp1.AND.ll_tmp2) THEN
432	zcpy(ji,jj) = 1.0_wp
433	wdvflt(ji,jj) = 1.0_wp
434	ELSE IF(ll_tmp3) THEN
435	! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here
436	zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) / &
437	& (sshn(ji,jj+1) - sshn(ji,jj)))
438	wdvflt(ji,jj) = 1.0_wp
439	ELSE
440	zcpy(ji,jj) = 0._wp
441	wdvflt(ji,jj) = 0.0_wp
442	END IF
443	END DO
444	END DO
445	CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp )
446	ENDIF
447
448
449	!jii=jidbg-nimpp+1;jjj=jjdbg-njmpp+1
450
451	! Surface value
452	DO jj = 2, jpjm1
453	DO ji = fs_2, fs_jpim1 ! vector opt.
454	! hydrostatic pressure gradient along s-surfaces
455	zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) &
456	& - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) )
457	zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) &
458	& - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) )
459	! s-coordinate pressure gradient correction
460	zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) &
461	& * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
462	zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) &
463	& * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
464
465
466	IF(ln_wd) THEN
467
468	zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
469	zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)
470	zuap = zuap * zcpx(ji,jj)
471	zvap = zvap * zcpy(ji,jj)
472	ENDIF
473
474	! add to the general momentum trend
475	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
476	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
477	END DO
478	END DO
479
480	! interior value (2=<jk=<jpkm1)
481	DO jk = 2, jpkm1
482	DO jj = 2, jpjm1
483	DO ji = fs_2, fs_jpim1 ! vector opt.
484	! hydrostatic pressure gradient along s-surfaces
485	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) &
486	& * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) &
487	& - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) )
488	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) &
489	& * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) &
490	& - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) )
491	! s-coordinate pressure gradient correction
492	zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) &
493	& * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj)
494	zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) &
495	& * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
496
497	IF(ln_wd) THEN
498	zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
499	zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)
500	zuap = zuap * zcpx(ji,jj)
501	zvap = zvap * zcpy(ji,jj)
502	ENDIF
503
504	! add to the general momentum trend
505	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
506	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
507	END DO
508	END DO
509	END DO
510	!
511	CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj )
512	IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
513	!
514	END SUBROUTINE hpg_sco
515
516	SUBROUTINE hpg_isf( kt )
517	!!---------------------------------------------------------------------
518	!! * ROUTINE hpg_sco *
519	!!
520	!! ** Method : s-coordinate case. Jacobian scheme.
521	!! The now hydrostatic pressure gradient at a given level, jk,
522	!! is computed by taking the vertical integral of the in-situ
523	!! density gradient along the model level from the suface to that
524	!! level. s-coordinates (ln_sco): a corrective term is added
525	!! to the horizontal pressure gradient :
526	!! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ]
527	!! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ]
528	!! add it to the general momentum trend (ua,va).
529	!! ua = ua - 1/e1u * zhpi
530	!! va = va - 1/e2v * zhpj
531	!! iceload is added and partial cell case are added to the top and bottom
532	!!
533	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
534	!!----------------------------------------------------------------------
535	INTEGER, INTENT(in) :: kt ! ocean time-step index
536	!!
537	INTEGER :: ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j ! dummy loop indices
538	REAL(wp) :: zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1 ! temporary scalars
539	REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj, zrhd
540	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop
541	REAL(wp), POINTER, DIMENSION(:,:) :: ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj
542	!!----------------------------------------------------------------------
543	!
544	CALL wrk_alloc( jpi,jpj, 2, ztstop)
545	CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd)
546	CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj)
547	!
548	IF( kt == nit000 ) THEN
549	IF(lwp) WRITE(numout,*)
550	IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf'
551	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used'
552	ENDIF
553
554	! Local constant initialization
555	zcoef0 = - grav * 0.5_wp
556	! To use density and not density anomaly
557	! IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume
558	! ELSE ; znad = 0._wp ! Fixed volume
559	! ENDIF
560	znad=1._wp
561	! iniitialised to 0. zhpi zhpi
562	zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp
563
564	!==================================================================================
565	!=====Compute iceload and contribution of the half first wet layer =================
566	!===================================================================================
567
568	! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)
569	ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp
570
571	! compute density of the water displaced by the ice shelf
572	zrhd = rhd ! save rhd
573	DO jk = 1, jpk
574	zdept(:,:)=gdept_1d(jk)
575	CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk))
576	END DO
577	WHERE ( tmask(:,:,:) == 1._wp)
578	rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd
579	END WHERE
580
581	! compute rhd at the ice/oce interface (ice shelf side)
582	CALL eos(ztstop,risfdep,zrhdtop_isf)
583
584	! compute rhd at the ice/oce interface (ocean side)
585	DO ji=1,jpi
586	DO jj=1,jpj
587	ikt=mikt(ji,jj)
588	ztstop(ji,jj,1)=tsn(ji,jj,ikt,1)
589	ztstop(ji,jj,2)=tsn(ji,jj,ikt,2)
590	END DO
591	END DO
592	CALL eos(ztstop,risfdep,zrhdtop_oce)
593	!
594	! Surface value + ice shelf gradient
595	! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v)
596	ziceload = 0._wp
597	DO jj = 1, jpj
598	DO ji = 1, jpi ! vector opt.
599	ikt=mikt(ji,jj)
600	ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1))
601	DO jk=2,ikt-1
602	ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) &
603	& * (1._wp - tmask(ji,jj,jk))
604	END DO
605	IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) &
606	& * ( risfdep(ji,jj) - gdept_1d(ikt-1) )
607	END DO
608	END DO
609	riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5
610	! compute zp from z=0 to first T wet point (correction due to zps not yet applied)
611	DO jj = 2, jpjm1
612	DO ji = fs_2, fs_jpim1 ! vector opt.
613	ikt=mikt(ji,jj) ; iktp1i=mikt(ji+1,jj); iktp1j=mikt(ji,jj+1)
614	! hydrostatic pressure gradient along s-surfaces and ice shelf pressure
615	! we assume ISF is in isostatic equilibrium
616	zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj ,iktp1i) &
617	& * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj ) ) &
618	& - 0.5_wp * fse3w(ji ,jj ,ikt ) &
619	& * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) &
620	& + ( ziceload(ji+1,jj) - ziceload(ji,jj)) )
621	zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji ,jj+1,iktp1j) &
622	& * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji ,jj+1) ) &
623	& - 0.5_wp * fse3w(ji ,jj ,ikt ) &
624	& * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) &
625	& + ( ziceload(ji,jj+1) - ziceload(ji,jj) ) )
626	! s-coordinate pressure gradient correction (=0 if z coordinate)
627	zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) &
628	& * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
629	zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) &
630	& * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
631	! add to the general momentum trend
632	ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1)
633	va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1)
634	END DO
635	END DO
636	!==================================================================================
637	!===== Compute partial cell contribution for the top cell =========================
638	!==================================================================================
639	DO jj = 2, jpjm1
640	DO ji = fs_2, fs_jpim1 ! vector opt.
641	iku = miku(ji,jj) ;
642	zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp
643	ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))
644	! u direction
645	IF ( iku .GT. 1 ) THEN
646	! case iku
647	zhpi(ji,jj,iku) = zcoef0 / e1u(ji,jj) * ze3wu &
648	& * ( rhd (ji+1,jj,iku) + rhd (ji,jj,iku) &
649	& + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )
650	! corrective term ( = 0 if z coordinate )
651	zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj)
652	! zhpi will be added in interior loop
653	ua(ji,jj,iku) = ua(ji,jj,iku) + zuap
654	! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure
655	IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj) = zhpi(ji,jj,iku)
656
657	! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)
658	zhpiint = zcoef0 / e1u(ji,jj) &
659	& * ( fse3w(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) &
660	& + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) &
661	& - fse3w(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) &
662	& + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) )
663	zhpi(ji,jj,iku+1) = zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint
664	END IF
665
666	! v direction
667	ikv = mikv(ji,jj)
668	ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))
669	IF ( ikv .GT. 1 ) THEN
670	! case ikv
671	zhpj(ji,jj,ikv) = zcoef0 / e2v(ji,jj) * ze3wv &
672	& * ( rhd(ji,jj+1,ikv) + rhd (ji,jj,ikv) &
673	& + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )
674	! corrective term ( = 0 if z coordinate )
675	zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj)
676	! zhpi will be added in interior loop
677	va(ji,jj,ikv) = va(ji,jj,ikv) + zvap
678	! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure
679	IF (mbkv(ji,jj) == ikv + 1) zpshpj(ji,jj) = zhpj(ji,jj,ikv)
680
681	! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)
682	zhpjint = zcoef0 / e2v(ji,jj) &
683	& * ( fse3w(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) &
684	& + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) &
685	& - fse3w(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) &
686	& + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) )
687	zhpj(ji,jj,ikv+1) = zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint
688	END IF
689	END DO
690	END DO
691
692	!==================================================================================
693	!===== Compute interior value =====================================================
694	!==================================================================================
695
696	DO jj = 2, jpjm1
697	DO ji = fs_2, fs_jpim1 ! vector opt.
698	iku=miku(ji,jj); ikv=mikv(ji,jj)
699	DO jk = 2, jpkm1
700	! hydrostatic pressure gradient along s-surfaces
701	! zhpi is masked for the first wet cell (contribution already done in the upper bloc)
702	zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1) &
703	& + zcoef0 / e1u(ji,jj) &
704	& * ( fse3w(ji+1,jj ,jk) * ( (rhd(ji+1,jj,jk ) + znad) &
705	& + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1) &
706	& - fse3w(ji ,jj ,jk) * ( (rhd(ji ,jj,jk ) + znad) &
707	& + (rhd(ji ,jj,jk-1) + znad) ) * tmask(ji ,jj,jk-1) )
708	! s-coordinate pressure gradient correction
709	! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)
710	zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) &
711	& * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1)
712	ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk)
713
714	! hydrostatic pressure gradient along s-surfaces
715	! zhpi is masked for the first wet cell (contribution already done in the upper bloc)
716	zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1) &
717	& + zcoef0 / e2v(ji,jj) &
718	& * ( fse3w(ji ,jj+1,jk) * ( (rhd(ji,jj+1,jk ) + znad) &
719	& + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1) &
720	& - fse3w(ji ,jj ,jk) * ( (rhd(ji,jj ,jk ) + znad) &
721	& + (rhd(ji,jj ,jk-1) + znad) ) * tmask(ji,jj ,jk-1) )
722	! s-coordinate pressure gradient correction
723	! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc)
724	zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) &
725	& * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1)
726	! add to the general momentum trend
727	va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk)
728	END DO
729	END DO
730	END DO
731
732	!==================================================================================
733	!===== Compute bottom cell contribution (partial cell) ============================
734	!==================================================================================
735
736	DO jj = 2, jpjm1
737	DO ji = 2, jpim1
738	iku = mbku(ji,jj)
739	ikv = mbkv(ji,jj)
740
741	IF (iku .GT. 1) THEN
742	! remove old value (interior case)
743	zuap = -zcoef0 * ( rhd (ji+1,jj ,iku) + rhd (ji,jj,iku) + 2._wp * znad ) &
744	& * ( fsde3w(ji+1,jj ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj)
745	ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap
746	! put new value
747	! -zpshpi to avoid double contribution of the partial step in the top layer
748	zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) / e1u(ji,jj)
749	zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj)
750	ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap
751	END IF
752	! v direction
753	IF (ikv .GT. 1) THEN
754	! remove old value (interior case)
755	zvap = -zcoef0 * ( rhd (ji ,jj+1,ikv) + rhd (ji,jj,ikv) + 2._wp * znad ) &
756	& * ( fsde3w(ji ,jj+1,ikv) - fsde3w(ji,jj,ikv) ) / e2v(ji,jj)
757	va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap
758	! put new value
759	! -zpshpj to avoid double contribution of the partial step in the top layer
760	zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) / e2v(ji,jj)
761	zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj)
762	va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap
763	END IF
764	END DO
765	END DO
766
767	! set back to original density value into the ice shelf cell (maybe useless because it is masked)
768	rhd = zrhd
769	!
770	CALL wrk_dealloc( jpi,jpj,2, ztstop)
771	CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd)
772	CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj)
773	!
774	END SUBROUTINE hpg_isf
775
776
777	SUBROUTINE hpg_djc( kt )
778	!!---------------------------------------------------------------------
779	!! * ROUTINE hpg_djc *
780	!!
781	!! ** Method : Density Jacobian with Cubic polynomial scheme
782	!!
783	!! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003
784	!!----------------------------------------------------------------------
785	INTEGER, INTENT(in) :: kt ! ocean time-step index
786	!!
787	INTEGER :: ji, jj, jk ! dummy loop indices
788	REAL(wp) :: zcoef0, zep, cffw ! temporary scalars
789	REAL(wp) :: z1_10, cffu, cffx ! " "
790	REAL(wp) :: z1_12, cffv, cffy ! " "
791	LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables
792	REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj
793	REAL(wp), POINTER, DIMENSION(:,:,:) :: dzx, dzy, dzz, dzu, dzv, dzw
794	REAL(wp), POINTER, DIMENSION(:,:,:) :: drhox, drhoy, drhoz, drhou, drhov, drhow
795	REAL(wp), POINTER, DIMENSION(:,:,:) :: rho_i, rho_j, rho_k
796	REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter
797	!!----------------------------------------------------------------------
798	!
799	CALL wrk_alloc( jpi, jpj, jpk, dzx , dzy , dzz , dzu , dzv , dzw )
800	CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow )
801	CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj )
802	IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
803	!
804	!!----------------------------------------------------------------------
805	!
806	!
807	IF(ln_wd) THEN
808	DO jj = 2, jpjm1
809	DO ji = 2, jpim1
810	ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) &
811	& .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) &
812	& > rn_wdmin1 + rn_wdmin2
813	ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +&
814	& rn_wdmin1 + rn_wdmin2
815
816	IF(ll_tmp1) THEN
817	zcpx(ji,jj) = 1.0_wp
818	ELSE IF(ll_tmp2) THEN
819	! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here
820	zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&
821	& (sshn(ji+1,jj) - sshn(ji,jj)))
822	ELSE
823	zcpx(ji,jj) = 0._wp
824	END IF
825
826	ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) &
827	& .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) &
828	& > rn_wdmin1 + rn_wdmin2
829	ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +&
830	& rn_wdmin1 + rn_wdmin2
831
832	IF(ll_tmp1) THEN
833	zcpy(ji,jj) = 1.0_wp
834	ELSE IF(ll_tmp2) THEN
835	! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here
836	zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /&
837	& (sshn(ji,jj+1) - sshn(ji,jj)))
838	ELSE
839	zcpy(ji,jj) = 0._wp
840	END IF
841	END DO
842	END DO
843	CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp )
844	ENDIF
845
846
847	IF( kt == nit000 ) THEN
848	IF(lwp) WRITE(numout,*)
849	IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend'
850	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme'
851	ENDIF
852
853	! Local constant initialization
854	zcoef0 = - grav * 0.5_wp
855	z1_10 = 1._wp / 10._wp
856	z1_12 = 1._wp / 12._wp
857
858	!----------------------------------------------------------------------------------------
859	! compute and store in provisional arrays elementary vertical and horizontal differences
860	!----------------------------------------------------------------------------------------
861
862	!!bug gm Not a true bug, but... dzz=e3w for dzx, dzy verify what it is really
863
864	DO jk = 2, jpkm1
865	DO jj = 2, jpjm1
866	DO ji = fs_2, fs_jpim1 ! vector opt.
867	drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1)
868	dzz (ji,jj,jk) = fsde3w(ji ,jj ,jk) - fsde3w(ji,jj,jk-1)
869	drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk )
870	dzx (ji,jj,jk) = fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk )
871	drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk )
872	dzy (ji,jj,jk) = fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk )
873	END DO
874	END DO
875	END DO
876
877	!-------------------------------------------------------------------------
878	! compute harmonic averages using eq. 5.18
879	!-------------------------------------------------------------------------
880	zep = 1.e-15
881
882	!!bug gm drhoz not defined at level 1 and used (jk-1 with jk=2)
883	!!bug gm idem for drhox, drhoy et ji=jpi and jj=jpj
884
885	DO jk = 2, jpkm1
886	DO jj = 2, jpjm1
887	DO ji = fs_2, fs_jpim1 ! vector opt.
888	cffw = 2._wp * drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1)
889
890	cffu = 2._wp * drhox(ji+1,jj ,jk) * drhox(ji,jj,jk )
891	cffx = 2._wp * dzx (ji+1,jj ,jk) * dzx (ji,jj,jk )
892
893	cffv = 2._wp * drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk )
894	cffy = 2._wp * dzy (ji ,jj+1,jk) * dzy (ji,jj,jk )
895
896	IF( cffw > zep) THEN
897	drhow(ji,jj,jk) = 2._wp * drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1) &
898	& / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) )
899	ELSE
900	drhow(ji,jj,jk) = 0._wp
901	ENDIF
902
903	dzw(ji,jj,jk) = 2._wp * dzz(ji,jj,jk) * dzz(ji,jj,jk-1) &
904	& / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) )
905
906	IF( cffu > zep ) THEN
907	drhou(ji,jj,jk) = 2._wp * drhox(ji+1,jj,jk) * drhox(ji,jj,jk) &
908	& / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) )
909	ELSE
910	drhou(ji,jj,jk ) = 0._wp
911	ENDIF
912
913	IF( cffx > zep ) THEN
914	dzu(ji,jj,jk) = 2._wp * dzx(ji+1,jj,jk) * dzx(ji,jj,jk) &
915	& / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) )
916	ELSE
917	dzu(ji,jj,jk) = 0._wp
918	ENDIF
919
920	IF( cffv > zep ) THEN
921	drhov(ji,jj,jk) = 2._wp * drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk) &
922	& / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) )
923	ELSE
924	drhov(ji,jj,jk) = 0._wp
925	ENDIF
926
927	IF( cffy > zep ) THEN
928	dzv(ji,jj,jk) = 2._wp * dzy(ji,jj+1,jk) * dzy(ji,jj,jk) &
929	& / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) )
930	ELSE
931	dzv(ji,jj,jk) = 0._wp
932	ENDIF
933
934	END DO
935	END DO
936	END DO
937
938	!----------------------------------------------------------------------------------
939	! apply boundary conditions at top and bottom using 5.36-5.37
940	!----------------------------------------------------------------------------------
941	drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:, 1 ) ) - 0.5_wp * drhow(:,:, 2 )
942	drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:, 1 ) ) - 0.5_wp * drhou(:,:, 2 )
943	drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:, 1 ) ) - 0.5_wp * drhov(:,:, 2 )
944
945	drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1)
946	drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1)
947	drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1)
948
949
950	!--------------------------------------------------------------
951	! Upper half of top-most grid box, compute and store
952	!-------------------------------------------------------------
953
954	!!bug gm : e3w-de3w = 0.5*e3w .... and de3w(2)-de3w(1)=e3w(2) .... to be verified
955	! true if de3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be
956
957	DO jj = 2, jpjm1
958	DO ji = fs_2, fs_jpim1 ! vector opt.
959	rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) ) &
960	& * ( rhd(ji,jj,1) &
961	& + 0.5_wp * ( rhd(ji,jj,2) - rhd(ji,jj,1) ) &
962	& * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) ) &
963	& / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) ) )
964	END DO
965	END DO
966
967	!!bug gm : here also, simplification is possible
968	!!bug gm : optimisation: 1/10 and 1/12 the division should be done before the loop
969
970	DO jk = 2, jpkm1
971	DO jj = 2, jpjm1
972	DO ji = fs_2, fs_jpim1 ! vector opt.
973
974	rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) &
975	& * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) ) &
976	& - grav * z1_10 * ( &
977	& ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) &
978	& * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) &
979	& - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) &
980	& * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) &
981	& )
982
983	rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) &
984	& * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) ) &
985	& - grav* z1_10 * ( &
986	& ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) &
987	& * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) &
988	& - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) &
989	& * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) &
990	& )
991
992	rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) &
993	& * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) ) &
994	& - grav* z1_10 * ( &
995	& ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) &
996	& * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) &
997	& - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) &
998	& * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) &
999	& )
1000
1001	END DO
1002	END DO
1003	END DO
1004	CALL lbc_lnk(rho_k,'W',1.)
1005	CALL lbc_lnk(rho_i,'U',1.)
1006	CALL lbc_lnk(rho_j,'V',1.)
1007
1008
1009	! ---------------
1010	! Surface value
1011	! ---------------
1012	DO jj = 2, jpjm1
1013	DO ji = fs_2, fs_jpim1 ! vector opt.
1014	zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj)
1015	zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj)
1016	IF(ln_wd) THEN
1017	zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj)
1018	zhpj(ji,jj,1) = zhpj(ji,jj,1) * zcpy(ji,jj)
1019	ENDIF
1020	! add to the general momentum trend
1021	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
1022	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
1023	END DO
1024	END DO
1025
1026	! ----------------
1027	! interior value (2=<jk=<jpkm1)
1028	! ----------------
1029	DO jk = 2, jpkm1
1030	DO jj = 2, jpjm1
1031	DO ji = fs_2, fs_jpim1 ! vector opt.
1032	! hydrostatic pressure gradient along s-surfaces
1033	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &
1034	& + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) ) &
1035	& - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) / e1u(ji,jj)
1036	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &
1037	& + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) &
1038	& -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) / e2v(ji,jj)
1039	IF(ln_wd) THEN
1040	zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
1041	zhpj(ji,jj,jk) = zhpj(ji,jj,jk) * zcpy(ji,jj)
1042	ENDIF
1043	! add to the general momentum trend
1044	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
1045	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
1046	END DO
1047	END DO
1048	END DO
1049	!
1050	CALL wrk_dealloc( jpi, jpj, jpk, dzx , dzy , dzz , dzu , dzv , dzw )
1051	CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow )
1052	CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj )
1053	IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
1054	!
1055	END SUBROUTINE hpg_djc
1056
1057
1058	SUBROUTINE hpg_prj( kt )
1059	!!---------------------------------------------------------------------
1060	!! * ROUTINE hpg_prj *
1061	!!
1062	!! ** Method : s-coordinate case.
1063	!! A Pressure-Jacobian horizontal pressure gradient method
1064	!! based on the constrained cubic-spline interpolation for
1065	!! all vertical coordinate systems
1066	!!
1067	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
1068	!!----------------------------------------------------------------------
1069	INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear
1070	INTEGER, INTENT(in) :: kt ! ocean time-step index
1071	!!
1072	INTEGER :: ji, jj, jk, jkk ! dummy loop indices
1073	REAL(wp) :: zcoef0, znad ! temporary scalars
1074	!!
1075	!! The local variables for the correction term
1076	INTEGER :: jk1, jis, jid, jjs, jjd
1077	LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables
1078	REAL(wp) :: zuijk, zvijk, zpwes, zpwed, zpnss, zpnsd, zdeps
1079	REAL(wp) :: zrhdt1
1080	REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2
1081	REAL(wp), POINTER, DIMENSION(:,:,:) :: zdept, zrhh
1082	REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp
1083	REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_n, zsshv_n
1084	REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy !W/D pressure filter
1085	!!----------------------------------------------------------------------
1086	!
1087	CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
1088	CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh )
1089	CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n )
1090	IF(ln_wd) CALL wrk_alloc( jpi,jpj, zcpx, zcpy )
1091	!
1092	IF( kt == nit000 ) THEN
1093	IF(lwp) WRITE(numout,*)
1094	IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend'
1095	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, cubic spline pressure Jacobian'
1096	ENDIF
1097
1098	!!----------------------------------------------------------------------
1099	! Local constant initialization
1100	zcoef0 = - grav
1101	znad = 0.0_wp
1102	IF( lk_vvl ) znad = 1._wp
1103
1104	IF(ln_wd) THEN
1105	DO jj = 2, jpjm1
1106	DO ji = 2, jpim1
1107	ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) &
1108	& .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) &
1109	& > rn_wdmin1 + rn_wdmin2
1110	ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) +&
1111	& rn_wdmin1 + rn_wdmin2
1112
1113	IF(ll_tmp1) THEN
1114	zcpx(ji,jj) = 1.0_wp
1115	ELSE IF(ll_tmp2) THEN
1116	! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here
1117	zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /&
1118	& (sshn(ji+1,jj) - sshn(ji,jj)))
1119	ELSE
1120	zcpx(ji,jj) = 0._wp
1121	END IF
1122
1123	ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) &
1124	& .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) &
1125	& > rn_wdmin1 + rn_wdmin2
1126	ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) +&
1127	& rn_wdmin1 + rn_wdmin2
1128
1129	IF(ll_tmp1.OR.ll_tmp2) THEN
1130	zcpy(ji,jj) = 1.0_wp
1131	ELSE IF(ll_tmp2) THEN
1132	! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here
1133	zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) /&
1134	& (sshn(ji,jj+1) - sshn(ji,jj)))
1135	ELSE
1136	zcpy(ji,jj) = 0._wp
1137	END IF
1138	END DO
1139	END DO
1140	CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp )
1141	ENDIF
1142
1143	! Clean 3-D work arrays
1144	zhpi(:,:,:) = 0._wp
1145	zrhh(:,:,:) = rhd(:,:,:)
1146
1147	! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate
1148	DO jj = 1, jpj
1149	DO ji = 1, jpi
1150	jk = mbathy(ji,jj)
1151	IF( jk <= 0 ) THEN; zrhh(ji,jj,:) = 0._wp
1152	ELSE IF(jk == 1) THEN; zrhh(ji,jj, jk+1:jpk) = rhd(ji,jj,jk)
1153	ELSE IF(jk < jpkm1) THEN
1154	DO jkk = jk+1, jpk
1155	zrhh(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1), &
1156	fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2))
1157	END DO
1158	ENDIF
1159	END DO
1160	END DO
1161
1162	! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)"
1163	DO jj = 1, jpj
1164	DO ji = 1, jpi
1165	zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad
1166	END DO
1167	END DO
1168
1169	DO jk = 2, jpk
1170	DO jj = 1, jpj
1171	DO ji = 1, jpi
1172	zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk)
1173	END DO
1174	END DO
1175	END DO
1176
1177	fsp(:,:,:) = zrhh (:,:,:)
1178	xsp(:,:,:) = zdept(:,:,:)
1179
1180	! Construct the vertical density profile with the
1181	! constrained cubic spline interpolation
1182	! rho(z) = asp + bspz + cspz^2 + dsp*z^3
1183	CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type)
1184
1185	! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)"
1186	DO jj = 2, jpj
1187	DO ji = 2, jpi
1188	zrhdt1 = zrhh(ji,jj,1) - interp3(zdept(ji,jj,1),asp(ji,jj,1), &
1189	bsp(ji,jj,1), csp(ji,jj,1), &
1190	dsp(ji,jj,1) ) * 0.25_wp * fse3w(ji,jj,1)
1191
1192	! assuming linear profile across the top half surface layer
1193	zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * zrhdt1
1194	END DO
1195	END DO
1196
1197	! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)"
1198	DO jk = 2, jpkm1
1199	DO jj = 2, jpj
1200	DO ji = 2, jpi
1201	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + &
1202	integ_spline(zdept(ji,jj,jk-1), zdept(ji,jj,jk),&
1203	asp(ji,jj,jk-1), bsp(ji,jj,jk-1), &
1204	csp(ji,jj,jk-1), dsp(ji,jj,jk-1))
1205	END DO
1206	END DO
1207	END DO
1208
1209	! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1)
1210
1211	! Prepare zsshu_n and zsshv_n
1212	DO jj = 2, jpjm1
1213	DO ji = 2, jpim1
1214	zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * &
1215	& r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp
1216	zsshv_n(ji,jj) = (e1e2v(ji,jj) * sshn(ji,jj) + e1e2v(ji+1, jj) * sshn(ji,jj+1)) * &
1217	& r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp
1218	END DO
1219	END DO
1220
1221	CALL lbc_lnk (zsshu_n, 'U', 1.)
1222	CALL lbc_lnk (zsshv_n, 'V', 1.)
1223
1224	DO jj = 2, jpjm1
1225	DO ji = 2, jpim1
1226	zu(ji,jj,1) = - ( fse3u(ji,jj,1) - zsshu_n(ji,jj) * znad)
1227	zv(ji,jj,1) = - ( fse3v(ji,jj,1) - zsshv_n(ji,jj) * znad)
1228	END DO
1229	END DO
1230
1231	DO jk = 2, jpkm1
1232	DO jj = 2, jpjm1
1233	DO ji = 2, jpim1
1234	zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u(ji,jj,jk)
1235	zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v(ji,jj,jk)
1236	END DO
1237	END DO
1238	END DO
1239
1240	DO jk = 1, jpkm1
1241	DO jj = 2, jpjm1
1242	DO ji = 2, jpim1
1243	zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk)
1244	zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk)
1245	END DO
1246	END DO
1247	END DO
1248
1249	DO jk = 1, jpkm1
1250	DO jj = 2, jpjm1
1251	DO ji = 2, jpim1
1252	zu(ji,jj,jk) = min(zu(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))
1253	zu(ji,jj,jk) = max(zu(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))
1254	zv(ji,jj,jk) = min(zv(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))
1255	zv(ji,jj,jk) = max(zv(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))
1256	END DO
1257	END DO
1258	END DO
1259
1260
1261	DO jk = 1, jpkm1
1262	DO jj = 2, jpjm1
1263	DO ji = 2, jpim1
1264	zpwes = 0._wp; zpwed = 0._wp
1265	zpnss = 0._wp; zpnsd = 0._wp
1266	zuijk = zu(ji,jj,jk)
1267	zvijk = zv(ji,jj,jk)
1268
1269	!!!!! for u equation
1270	IF( jk <= mbku(ji,jj) ) THEN
1271	IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN
1272	jis = ji + 1; jid = ji
1273	ELSE
1274	jis = ji; jid = ji +1
1275	ENDIF
1276
1277	! integrate the pressure on the shallow side
1278	jk1 = jk
1279	DO WHILE ( -zdept(jis,jj,jk1) > zuijk )
1280	IF( jk1 == mbku(ji,jj) ) THEN
1281	zuijk = -zdept(jis,jj,jk1)
1282	EXIT
1283	ENDIF
1284	zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk)
1285	zpwes = zpwes + &
1286	integ_spline(zdept(jis,jj,jk1), zdeps, &
1287	asp(jis,jj,jk1), bsp(jis,jj,jk1), &
1288	csp(jis,jj,jk1), dsp(jis,jj,jk1))
1289	jk1 = jk1 + 1
1290	END DO
1291
1292	! integrate the pressure on the deep side
1293	jk1 = jk
1294	DO WHILE ( -zdept(jid,jj,jk1) < zuijk )
1295	IF( jk1 == 1 ) THEN
1296	zdeps = zdept(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad)
1297	zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), &
1298	bsp(jid,jj,1), csp(jid,jj,1), &
1299	dsp(jid,jj,1)) * zdeps
1300	zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps
1301	EXIT
1302	ENDIF
1303	zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk)
1304	zpwed = zpwed + &
1305	integ_spline(zdeps, zdept(jid,jj,jk1), &
1306	asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), &
1307	csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) )
1308	jk1 = jk1 - 1
1309	END DO
1310
1311	! update the momentum trends in u direction
1312
1313	zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))
1314	IF( lk_vvl ) THEN
1315	zdpdx2 = zcoef0 / e1u(ji,jj) * &
1316	( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) )
1317	ELSE
1318	zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed)
1319	ENDIF
1320
1321	IF(ln_wd) THEN
1322	zdpdx1 = zdpdx1 * zcpx(ji,jj)
1323	zdpdx2 = zdpdx2 * zcpx(ji,jj)
1324	ENDIF
1325	ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk)
1326	ENDIF
1327
1328	!!!!! for v equation
1329	IF( jk <= mbkv(ji,jj) ) THEN
1330	IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN
1331	jjs = jj + 1; jjd = jj
1332	ELSE
1333	jjs = jj ; jjd = jj + 1
1334	ENDIF
1335
1336	! integrate the pressure on the shallow side
1337	jk1 = jk
1338	DO WHILE ( -zdept(ji,jjs,jk1) > zvijk )
1339	IF( jk1 == mbkv(ji,jj) ) THEN
1340	zvijk = -zdept(ji,jjs,jk1)
1341	EXIT
1342	ENDIF
1343	zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk)
1344	zpnss = zpnss + &
1345	integ_spline(zdept(ji,jjs,jk1), zdeps, &
1346	asp(ji,jjs,jk1), bsp(ji,jjs,jk1), &
1347	csp(ji,jjs,jk1), dsp(ji,jjs,jk1) )
1348	jk1 = jk1 + 1
1349	END DO
1350
1351	! integrate the pressure on the deep side
1352	jk1 = jk
1353	DO WHILE ( -zdept(ji,jjd,jk1) < zvijk )
1354	IF( jk1 == 1 ) THEN
1355	zdeps = zdept(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad)
1356	zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), &
1357	bsp(ji,jjd,1), csp(ji,jjd,1), &
1358	dsp(ji,jjd,1) ) * zdeps
1359	zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps
1360	EXIT
1361	ENDIF
1362	zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk)
1363	zpnsd = zpnsd + &
1364	integ_spline(zdeps, zdept(ji,jjd,jk1), &
1365	asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), &
1366	csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
1367	jk1 = jk1 - 1
1368	END DO
1369
1370
1371	! update the momentum trends in v direction
1372
1373	zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))
1374	IF( lk_vvl ) THEN
1375	zdpdy2 = zcoef0 / e2v(ji,jj) * &
1376	( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) )
1377	ELSE
1378	zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )
1379	ENDIF
1380
1381	IF(ln_wd) THEN
1382	zdpdy1 = zdpdy1 * zcpy(ji,jj)
1383	zdpdy2 = zdpdy2 * zcpy(ji,jj)
1384	ENDIF
1385
1386	va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2) * vmask(ji,jj,jk)
1387	ENDIF
1388
1389
1390	END DO
1391	END DO
1392	END DO
1393	!
1394	CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )
1395	CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh )
1396	CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n )
1397	IF(ln_wd) CALL wrk_dealloc( jpi,jpj, zcpx, zcpy )
1398	!
1399	END SUBROUTINE hpg_prj
1400
1401
1402	SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type)
1403	!!----------------------------------------------------------------------
1404	!! * ROUTINE cspline *
1405	!!
1406	!! ** Purpose : constrained cubic spline interpolation
1407	!!
1408	!! ** Method : f(x) = asp + bspx + cspx^2 + dsp*x^3
1409	!!
1410	!! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation
1411	!!----------------------------------------------------------------------
1412	IMPLICIT NONE
1413	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: fsp, xsp ! value and coordinate
1414	REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of
1415	! the interpoated function
1416	INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline
1417	! 2: Linear
1418	!
1419	INTEGER :: ji, jj, jk ! dummy loop indices
1420	INTEGER :: jpi, jpj, jpkm1
1421	REAL(wp) :: zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp
1422	REAL(wp) :: zdxtmp1, zdxtmp2, zalpha
1423	REAL(wp) :: zdf(size(fsp,3))
1424	!!----------------------------------------------------------------------
1425
1426	jpi = size(fsp,1)
1427	jpj = size(fsp,2)
1428	jpkm1 = size(fsp,3) - 1
1429
1430
1431	IF (polynomial_type == 1) THEN ! Constrained Cubic Spline
1432	DO ji = 1, jpi
1433	DO jj = 1, jpj
1434	!!Fritsch&Butland's method, 1984 (preferred, but more computation)
1435	! DO jk = 2, jpkm1-1
1436	! zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1)
1437	! zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1438	! zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1
1439	! zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2
1440	!
1441	! zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp
1442	!
1443	! IF(zdf1 * zdf2 <= 0._wp) THEN
1444	! zdf(jk) = 0._wp
1445	! ELSE
1446	! zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 )
1447	! ENDIF
1448	! END DO
1449
1450	!!Simply geometric average
1451	DO jk = 2, jpkm1-1
1452	zdf1 = (fsp(ji,jj,jk) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk) - xsp(ji,jj,jk-1))
1453	zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk))
1454
1455	IF(zdf1 * zdf2 <= 0._wp) THEN
1456	zdf(jk) = 0._wp
1457	ELSE
1458	zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2)
1459	ENDIF
1460	END DO
1461
1462	zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / &
1463	& ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2)
1464	zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / &
1465	& ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - &
1466	& 0.5_wp * zdf(jpkm1 - 1)
1467
1468	DO jk = 1, jpkm1 - 1
1469	zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1470	ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp
1471	ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp
1472	zddf1 = -2._wp * ztmp1 + ztmp2
1473	ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp
1474	zddf2 = 2._wp * ztmp1 - ztmp2
1475
1476	dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp
1477	csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp
1478	bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - &
1479	& csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - &
1480	& dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - &
1481	& xsp(ji,jj,jk+1) * xsp(ji,jj,jk))
1482	asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + &
1483	& (xsp(ji,jj,jk) * (csp(ji,jj,jk) + &
1484	& dsp(ji,jj,jk) * xsp(ji,jj,jk))))
1485	END DO
1486	END DO
1487	END DO
1488
1489	ELSE IF (polynomial_type == 2) THEN ! Linear
1490	DO ji = 1, jpi
1491	DO jj = 1, jpj
1492	DO jk = 1, jpkm1-1
1493	zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1494	ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk)
1495
1496	dsp(ji,jj,jk) = 0._wp
1497	csp(ji,jj,jk) = 0._wp
1498	bsp(ji,jj,jk) = ztmp1 / zdxtmp
1499	asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk)
1500	END DO
1501	END DO
1502	END DO
1503
1504	ELSE
1505	CALL ctl_stop( 'invalid polynomial type in cspline' )
1506	ENDIF
1507
1508	END SUBROUTINE cspline
1509
1510
1511	FUNCTION interp1(x, xl, xr, fl, fr) RESULT(f)
1512	!!----------------------------------------------------------------------
1513	!! * ROUTINE interp1 *
1514	!!
1515	!! ** Purpose : 1-d linear interpolation
1516	!!
1517	!! ** Method : interpolation is straight forward
1518	!! extrapolation is also permitted (no value limit)
1519	!!----------------------------------------------------------------------
1520	IMPLICIT NONE
1521	REAL(wp), INTENT(in) :: x, xl, xr, fl, fr
1522	REAL(wp) :: f ! result of the interpolation (extrapolation)
1523	REAL(wp) :: zdeltx
1524	!!----------------------------------------------------------------------
1525
1526	zdeltx = xr - xl
1527	IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN
1528	f = 0.5_wp * (fl + fr)
1529	ELSE
1530	f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx
1531	ENDIF
1532
1533	END FUNCTION interp1
1534
1535
1536	FUNCTION interp2(x, a, b, c, d) RESULT(f)
1537	!!----------------------------------------------------------------------
1538	!! * ROUTINE interp1 *
1539	!!
1540	!! ** Purpose : 1-d constrained cubic spline interpolation
1541	!!
1542	!! ** Method : cubic spline interpolation
1543	!!
1544	!!----------------------------------------------------------------------
1545	IMPLICIT NONE
1546	REAL(wp), INTENT(in) :: x, a, b, c, d
1547	REAL(wp) :: f ! value from the interpolation
1548	!!----------------------------------------------------------------------
1549
1550	f = a + x* ( b + x * ( c + d * x ) )
1551
1552	END FUNCTION interp2
1553
1554
1555	FUNCTION interp3(x, a, b, c, d) RESULT(f)
1556	!!----------------------------------------------------------------------
1557	!! * ROUTINE interp1 *
1558	!!
1559	!! ** Purpose : Calculate the first order of deriavtive of
1560	!! a cubic spline function y=a+bx+cx^2+d*x^3
1561	!!
1562	!! ** Method : f=dy/dx=b+2cx+3dx^2
1563	!!
1564	!!----------------------------------------------------------------------
1565	IMPLICIT NONE
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	IMPLICIT NONE
1585	REAL(wp), INTENT(in) :: xl, xr, a, b, c, d
1586	REAL(wp) :: za1, za2, za3
1587	REAL(wp) :: f ! integration result
1588	!!----------------------------------------------------------------------
1589
1590	za1 = 0.5_wp * b
1591	za2 = c / 3.0_wp
1592	za3 = 0.25_wp * d
1593
1594	f = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - &
1595	& xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) )
1596
1597	END FUNCTION integ_spline
1598
1599	!!======================================================================
1600	END MODULE dynhpg
1601

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

Download in other formats: