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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

dynhpg.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynhpg.F90 @ 10946

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

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

Download in other formats: