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/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 4747

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

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

Download in other formats: