Context Navigation

dynhpg.F90 @ 8403

Last change on this file since 8403 was 8403, checked in by deazer, 7 years ago

Add in ROMS WAD option ln_rwd+changes for implicit Bed Friction for ln_wd option
Note no ramp placed on ROMS bed friction yet
CS15mini case added as a Test CASE
at this revision AMM15 with Pure sigma coords barotorpic runs for 4 days without failure
in with ROMS option with 20cm min deoth and 50 vertical levels
Both run for CS15mini
In real domains nothing done on reference level yet so real domains
must have not negative depth points yet.
But a basic test has been done in WAD channel test cases (WAD7)

No changes in Main line source yet. See the MY_SRC sub dir of CS15 and TEST_CASES/WAD
for actual code changes.

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

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

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

Context Navigation

source: branches/UKMO/ROMS_WAD_7832/NEMOGCM/CONFIG/CS15mini/MY_SRC/dynhpg.F90 @ 8403

Download in other formats: