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

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

dynhpg.F90 in branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 3074

Last change on this file since 3074 was 3074, checked in by acc, 12 years ago
Branch dev_NOC_2011_MERGE. #874. Step 11. Merge in changes from the dev_r2802_NOCL_prjhpg. Partially imposed coding standards (more could be done); added missing namelist changes and added first attempt at documentation
Property svn:keywords set to `Id`
File size: 74.2 KB

Line
1	MODULE dynhpg
2	!!======================================================================
3	!! * MODULE dynhpg *
4	!! Ocean dynamics: hydrostatic pressure gradient trend
5	!!======================================================================
6	!! History : OPA ! 1987-09 (P. Andrich, M.-A. Foujols) hpg_zco: Original code
7	!! 5.0 ! 1991-11 (G. Madec)
8	!! 7.0 ! 1996-01 (G. Madec) hpg_sco: Original code for s-coordinates
9	!! 8.0 ! 1997-05 (G. Madec) split dynber into dynkeg and dynhpg
10	!! 8.5 ! 2002-07 (G. Madec) F90: Free form and module
11	!! 8.5 ! 2002-08 (A. Bozec) hpg_zps: Original code
12	!! NEMO 1.0 ! 2005-10 (A. Beckmann, B.W. An) various s-coordinate options
13	!! ! Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot
14	!! - ! 2005-11 (G. Madec) style & small optimisation
15	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
16	!!----------------------------------------------------------------------
17
18	!!----------------------------------------------------------------------
19	!! dyn_hpg : update the momentum trend with the now horizontal
20	!! gradient of the hydrostatic pressure
21	!! dyn_hpg_init : initialisation and control of options
22	!! hpg_zco : z-coordinate scheme
23	!! hpg_zps : z-coordinate plus partial steps (interpolation)
24	!! hpg_sco : s-coordinate (standard jacobian formulation)
25	!! hpg_hel : s-coordinate (helsinki modification)
26	!! hpg_wdj : s-coordinate (weighted density jacobian)
27	!! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial)
28	!! hpg_rot : s-coordinate (ROTated axes scheme)
29	!! hpg_prj : s-coordinate (Pressure Jacobian with Cubic polynomial)
30	!!----------------------------------------------------------------------
31	USE oce ! ocean dynamics and tracers
32	USE dom_oce ! ocean space and time domain
33	USE phycst ! physical constants
34	USE trdmod ! ocean dynamics trends
35	USE trdmod_oce ! ocean variables trends
36	USE in_out_manager ! I/O manager
37	USE prtctl ! Print control
38	USE lbclnk ! lateral boundary condition
39	USE lib_mpp ! MPP library
40
41	IMPLICIT NONE
42	PRIVATE
43
44	PUBLIC dyn_hpg ! routine called by step module
45	PUBLIC dyn_hpg_init ! routine called by opa module
46
47	! !!* Namelist namdyn_hpg : hydrostatic pressure gradient
48	LOGICAL , PUBLIC :: ln_hpg_zco = .TRUE. !: z-coordinate - full steps
49	LOGICAL , PUBLIC :: ln_hpg_zps = .FALSE. !: z-coordinate - partial steps (interpolation)
50	LOGICAL , PUBLIC :: ln_hpg_sco = .FALSE. !: s-coordinate (standard jacobian formulation)
51	LOGICAL , PUBLIC :: ln_hpg_hel = .FALSE. !: s-coordinate (helsinki modification)
52	LOGICAL , PUBLIC :: ln_hpg_wdj = .FALSE. !: s-coordinate (weighted density jacobian)
53	LOGICAL , PUBLIC :: ln_hpg_djc = .FALSE. !: s-coordinate (Density Jacobian with Cubic polynomial)
54	LOGICAL , PUBLIC :: ln_hpg_rot = .FALSE. !: s-coordinate (ROTated axes scheme)
55	LOGICAL , PUBLIC :: ln_hpg_prj = .FALSE. !: s-coordinate (Pressure Jacobian scheme)
56	REAL(wp), PUBLIC :: rn_gamma = 0._wp !: weighting coefficient
57	LOGICAL , PUBLIC :: ln_dynhpg_imp = .FALSE. !: semi-implicite hpg flag
58
59	INTEGER :: nhpg = 0 ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags)
60
61	!! * Substitutions
62	# include "domzgr_substitute.h90"
63	# include "vectopt_loop_substitute.h90"
64	!!----------------------------------------------------------------------
65	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
66	!! $Id$
67	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
68	!!----------------------------------------------------------------------
69	CONTAINS
70
71	SUBROUTINE dyn_hpg( kt )
72	!!---------------------------------------------------------------------
73	!! * ROUTINE dyn_hpg *
74	!!
75	!! ** Method : Call the hydrostatic pressure gradient routine
76	!! using the scheme defined in the namelist
77	!!
78	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
79	!! - Save the trend (l_trddyn=T)
80	!!----------------------------------------------------------------------
81	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
82	USE wrk_nemo, ONLY: ztrdu => wrk_3d_1 , ztrdv => wrk_3d_2 ! 3D workspace
83	!!
84	INTEGER, INTENT(in) :: kt ! ocean time-step index
85	!!----------------------------------------------------------------------
86	!
87	IF( wrk_in_use(3, 1,2) ) THEN
88	CALL ctl_stop('dyn_hpg: requested workspace arrays are unavailable') ; RETURN
89	ENDIF
90	!
91	IF( l_trddyn ) THEN ! Temporary saving of ua and va trends (l_trddyn)
92	ztrdu(:,:,:) = ua(:,:,:)
93	ztrdv(:,:,:) = va(:,:,:)
94	ENDIF
95	!
96	SELECT CASE ( nhpg ) ! Hydrastatic pressure gradient computation
97	CASE ( 0 ) ; CALL hpg_zco ( kt ) ! z-coordinate
98	CASE ( 1 ) ; CALL hpg_zps ( kt ) ! z-coordinate plus partial steps (interpolation)
99	CASE ( 2 ) ; CALL hpg_sco ( kt ) ! s-coordinate (standard jacobian formulation)
100	CASE ( 3 ) ; CALL hpg_hel ( kt ) ! s-coordinate (helsinki modification)
101	CASE ( 4 ) ; CALL hpg_wdj ( kt ) ! s-coordinate (weighted density jacobian)
102	CASE ( 5 ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial)
103	CASE ( 6 ) ; CALL hpg_rot ( kt ) ! s-coordinate (ROTated axes scheme)
104	CASE ( 7 ) ; CALL hpg_prj ( kt ) ! s-coordinate (Pressure Jacobian scheme)
105	END SELECT
106	!
107	IF( l_trddyn ) THEN ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
108	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
109	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
110	CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )
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( wrk_not_released(3, 1,2) ) CALL ctl_stop('dyn_hpg: failed to release workspace arrays')
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	!!
133	NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_hel, &
134	& ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, ln_hpg_prj, &
135	& rn_gamma , ln_dynhpg_imp
136	!!----------------------------------------------------------------------
137	!
138	REWIND( numnam ) ! Read Namelist namdyn_hpg
139	READ ( numnam, namdyn_hpg )
140	!
141	IF(lwp) THEN ! Control print
142	WRITE(numout,*)
143	WRITE(numout,*) 'dyn_hpg_init : hydrostatic pressure gradient initialisation'
144	WRITE(numout,*) '~~~~~~~~~~~~'
145	WRITE(numout,*) ' Namelist namdyn_hpg : choice of hpg scheme'
146	WRITE(numout,*) ' z-coord. - full steps ln_hpg_zco = ', ln_hpg_zco
147	WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps
148	WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco
149	WRITE(numout,*) ' s-coord. (helsinki modification) ln_hpg_hel = ', ln_hpg_hel
150	WRITE(numout,*) ' s-coord. (weighted density jacobian) ln_hpg_wdj = ', ln_hpg_wdj
151	WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc
152	WRITE(numout,*) ' s-coord. (ROTated axes scheme) ln_hpg_rot = ', ln_hpg_rot
153	WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj
154	WRITE(numout,*) ' weighting coeff. (wdj scheme) rn_gamma = ', rn_gamma
155	WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp
156	ENDIF
157	!
158	IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) ) &
159	& CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:&
160	& the standard jacobian formulation hpg_sco or &
161	& the pressure jacobian formulation hpg_prj')
162	!
163	! ! Set nhpg from ln_hpg_... flags
164	IF( ln_hpg_zco ) nhpg = 0
165	IF( ln_hpg_zps ) nhpg = 1
166	IF( ln_hpg_sco ) nhpg = 2
167	IF( ln_hpg_hel ) nhpg = 3
168	IF( ln_hpg_wdj ) nhpg = 4
169	IF( ln_hpg_djc ) nhpg = 5
170	IF( ln_hpg_rot ) nhpg = 6
171	IF( ln_hpg_prj ) nhpg = 7
172	!
173	! ! Consitency check
174	ioptio = 0
175	IF( ln_hpg_zco ) ioptio = ioptio + 1
176	IF( ln_hpg_zps ) ioptio = ioptio + 1
177	IF( ln_hpg_sco ) ioptio = ioptio + 1
178	IF( ln_hpg_hel ) ioptio = ioptio + 1
179	IF( ln_hpg_wdj ) ioptio = ioptio + 1
180	IF( ln_hpg_djc ) ioptio = ioptio + 1
181	IF( ln_hpg_rot ) ioptio = ioptio + 1
182	IF( ln_hpg_prj ) ioptio = ioptio + 1
183	IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
184	!
185	END SUBROUTINE dyn_hpg_init
186
187
188	SUBROUTINE hpg_zco( kt )
189	!!---------------------------------------------------------------------
190	!! * ROUTINE hpg_zco *
191	!!
192	!! ** Method : z-coordinate case, levels are horizontal surfaces.
193	!! The now hydrostatic pressure gradient at a given level, jk,
194	!! is computed by taking the vertical integral of the in-situ
195	!! density gradient along the model level from the suface to that
196	!! level: zhpi = grav .....
197	!! zhpj = grav .....
198	!! add it to the general momentum trend (ua,va).
199	!! ua = ua - 1/e1u * zhpi
200	!! va = va - 1/e2v * zhpj
201	!!
202	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
203	!!----------------------------------------------------------------------
204	USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
205	!!
206	INTEGER, INTENT(in) :: kt ! ocean time-step index
207	!!
208	INTEGER :: ji, jj, jk ! dummy loop indices
209	REAL(wp) :: zcoef0, zcoef1 ! temporary scalars
210	!!----------------------------------------------------------------------
211
212	IF( kt == nit000 ) THEN
213	IF(lwp) WRITE(numout,*)
214	IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend'
215	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case '
216	ENDIF
217
218	zcoef0 = - grav * 0.5_wp ! Local constant initialization
219
220	! Surface value
221	DO jj = 2, jpjm1
222	DO ji = fs_2, fs_jpim1 ! vector opt.
223	zcoef1 = zcoef0 * fse3w(ji,jj,1)
224	! hydrostatic pressure gradient
225	zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
226	zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
227	! add to the general momentum trend
228	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
229	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
230	END DO
231	END DO
232	!
233	! interior value (2=<jk=<jpkm1)
234	DO jk = 2, jpkm1
235	DO jj = 2, jpjm1
236	DO ji = fs_2, fs_jpim1 ! vector opt.
237	zcoef1 = zcoef0 * fse3w(ji,jj,jk)
238	! hydrostatic pressure gradient
239	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &
240	& + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) &
241	& - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) / e1u(ji,jj)
242
243	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &
244	& + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) &
245	& - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) / e2v(ji,jj)
246	! add to the general momentum trend
247	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
248	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
249	END DO
250	END DO
251	END DO
252	!
253	END SUBROUTINE hpg_zco
254
255
256	SUBROUTINE hpg_zps( kt )
257	!!---------------------------------------------------------------------
258	!! * ROUTINE hpg_zps *
259	!!
260	!! ** Method : z-coordinate plus partial steps case. blahblah...
261	!!
262	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
263	!!----------------------------------------------------------------------
264	USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
265	!!
266	INTEGER, INTENT(in) :: kt ! ocean time-step index
267	!!
268	INTEGER :: ji, jj, jk ! dummy loop indices
269	INTEGER :: iku, ikv ! temporary integers
270	REAL(wp) :: zcoef0, zcoef1, zcoef2, zcoef3 ! temporary scalars
271	!!----------------------------------------------------------------------
272
273	IF( kt == nit000 ) THEN
274	IF(lwp) WRITE(numout,*)
275	IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend'
276	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate with partial steps - vector optimization'
277	ENDIF
278
279	! Local constant initialization
280	zcoef0 = - grav * 0.5_wp
281
282	! Surface value (also valid in partial step case)
283	DO jj = 2, jpjm1
284	DO ji = fs_2, fs_jpim1 ! vector opt.
285	zcoef1 = zcoef0 * fse3w(ji,jj,1)
286	! hydrostatic pressure gradient
287	zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
288	zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
289	! add to the general momentum trend
290	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
291	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
292	END DO
293	END DO
294
295	! interior value (2=<jk=<jpkm1)
296	DO jk = 2, jpkm1
297	DO jj = 2, jpjm1
298	DO ji = fs_2, fs_jpim1 ! vector opt.
299	zcoef1 = zcoef0 * fse3w(ji,jj,jk)
300	! hydrostatic pressure gradient
301	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &
302	& + zcoef1 * ( ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) &
303	& - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) / e1u(ji,jj)
304
305	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &
306	& + zcoef1 * ( ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) &
307	& - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) / e2v(ji,jj)
308	! add to the general momentum trend
309	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
310	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
311	END DO
312	END DO
313	END DO
314
315	! partial steps correction at the last level (use gru & grv computed in zpshde.F90)
316	# if defined key_vectopt_loop
317	jj = 1
318	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
319	# else
320	DO jj = 2, jpjm1
321	DO ji = 2, jpim1
322	# endif
323	iku = mbku(ji,jj)
324	ikv = mbkv(ji,jj)
325	zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj ,iku) )
326	zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji ,jj+1,ikv) )
327	IF( iku > 1 ) THEN ! on i-direction (level 2 or more)
328	ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value
329	zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one
330	& + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj)
331	ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend
332	ENDIF
333	IF( ikv > 1 ) THEN ! on j-direction (level 2 or more)
334	va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value
335	zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one
336	& + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj)
337	va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend
338	ENDIF
339	# if ! defined key_vectopt_loop
340	END DO
341	# endif
342	END DO
343	!
344	END SUBROUTINE hpg_zps
345
346
347	SUBROUTINE hpg_sco( kt )
348	!!---------------------------------------------------------------------
349	!! * ROUTINE hpg_sco *
350	!!
351	!! ** Method : s-coordinate case. Jacobian scheme.
352	!! The now hydrostatic pressure gradient at a given level, jk,
353	!! is computed by taking the vertical integral of the in-situ
354	!! density gradient along the model level from the suface to that
355	!! level. s-coordinates (ln_sco): a corrective term is added
356	!! to the horizontal pressure gradient :
357	!! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ]
358	!! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ]
359	!! add it to the general momentum trend (ua,va).
360	!! ua = ua - 1/e1u * zhpi
361	!! va = va - 1/e2v * zhpj
362	!!
363	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
364	!!----------------------------------------------------------------------
365	USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
366	!!
367	INTEGER, INTENT(in) :: kt ! ocean time-step index
368	!!
369	INTEGER :: ji, jj, jk ! dummy loop indices
370	REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars
371	!!----------------------------------------------------------------------
372
373	IF( kt == nit000 ) THEN
374	IF(lwp) WRITE(numout,*)
375	IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
376	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used'
377	ENDIF
378
379	! Local constant initialization
380	zcoef0 = - grav * 0.5_wp
381	! To use density and not density anomaly
382	IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume
383	ELSE ; znad = 0._wp ! Fixed volume
384	ENDIF
385
386	! Surface value
387	DO jj = 2, jpjm1
388	DO ji = fs_2, fs_jpim1 ! vector opt.
389	! hydrostatic pressure gradient along s-surfaces
390	zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) &
391	& - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) )
392	zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) &
393	& - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) )
394	! s-coordinate pressure gradient correction
395	zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) &
396	& * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
397	zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) &
398	& * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
399	! add to the general momentum trend
400	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
401	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
402	END DO
403	END DO
404
405	! interior value (2=<jk=<jpkm1)
406	DO jk = 2, jpkm1
407	DO jj = 2, jpjm1
408	DO ji = fs_2, fs_jpim1 ! vector opt.
409	! hydrostatic pressure gradient along s-surfaces
410	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) &
411	& * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) &
412	& - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) )
413	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) &
414	& * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) &
415	& - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) )
416	! s-coordinate pressure gradient correction
417	zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) &
418	& * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj)
419	zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) &
420	& * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
421	! add to the general momentum trend
422	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
423	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
424	END DO
425	END DO
426	END DO
427	!
428	END SUBROUTINE hpg_sco
429
430
431	SUBROUTINE hpg_hel( kt )
432	!!---------------------------------------------------------------------
433	!! * ROUTINE hpg_hel *
434	!!
435	!! ** Method : s-coordinate case.
436	!! The now hydrostatic pressure gradient at a given level
437	!! jk is computed by taking the vertical integral of the in-situ
438	!! density gradient along the model level from the suface to that
439	!! level. s-coordinates (ln_sco): a corrective term is added
440	!! to the horizontal pressure gradient :
441	!! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ]
442	!! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ]
443	!! add it to the general momentum trend (ua,va).
444	!! ua = ua - 1/e1u * zhpi
445	!! va = va - 1/e2v * zhpj
446	!!
447	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
448	!! - Save the trend (l_trddyn=T)
449	!!----------------------------------------------------------------------
450	USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
451	!!
452	INTEGER, INTENT(in) :: kt ! ocean time-step index
453	!!
454	INTEGER :: ji, jj, jk ! dummy loop indices
455	REAL(wp) :: zcoef0, zuap, zvap ! temporary scalars
456	!!----------------------------------------------------------------------
457
458	IF( kt == nit000 ) THEN
459	IF(lwp) WRITE(numout,*)
460	IF(lwp) WRITE(numout,*) 'dyn:hpg_hel : hydrostatic pressure gradient trend'
461	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, helsinki modified scheme'
462	ENDIF
463
464	! Local constant initialization
465	zcoef0 = - grav * 0.5_wp
466
467	! Surface value
468	DO jj = 2, jpjm1
469	DO ji = fs_2, fs_jpim1 ! vector opt.
470	! hydrostatic pressure gradient along s-surfaces
471	zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj ,1) * rhd(ji+1,jj ,1) &
472	& - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) )
473	zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji ,jj+1,1) * rhd(ji ,jj+1,1) &
474	& - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) )
475	! s-coordinate pressure gradient correction
476	zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) &
477	& * ( fsdept(ji+1,jj,1) - fsdept(ji,jj,1) ) / e1u(ji,jj)
478	zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) &
479	& * ( fsdept(ji,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj)
480	! add to the general momentum trend
481	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
482	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
483	END DO
484	END DO
485	!
486	! interior value (2=<jk=<jpkm1)
487	DO jk = 2, jpkm1
488	DO jj = 2, jpjm1
489	DO ji = fs_2, fs_jpim1 ! vector opt.
490	! hydrostatic pressure gradient along s-surfaces
491	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &
492	& + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk ) * rhd(ji+1,jj,jk) &
493	& -fse3t(ji ,jj,jk ) * rhd(ji ,jj,jk) ) &
494	& + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1) &
495	& -fse3t(ji ,jj,jk-1) * rhd(ji ,jj,jk-1) )
496	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &
497	& + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk) &
498	& -fse3t(ji,jj ,jk ) * rhd(ji,jj, jk) ) &
499	& + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) &
500	& -fse3t(ji,jj ,jk-1) * rhd(ji,jj, jk-1) )
501	! s-coordinate pressure gradient correction
502	zuap = - zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) &
503	& * ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)
504	zvap = - zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) &
505	& * ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)
506	! add to the general momentum trend
507	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
508	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
509	END DO
510	END DO
511	END DO
512	!
513	END SUBROUTINE hpg_hel
514
515
516	SUBROUTINE hpg_wdj( kt )
517	!!---------------------------------------------------------------------
518	!! * ROUTINE hpg_wdj *
519	!!
520	!! ** Method : Weighted Density Jacobian (wdj) scheme (song 1998)
521	!! The weighting coefficients from the namelist parameter rn_gamma
522	!! (alpha=0.5-rn_gamma ; beta=1-alpha=0.5+rn_gamma
523	!!
524	!! Reference : Song, Mon. Wea. Rev., 126, 3213-3230, 1998.
525	!!----------------------------------------------------------------------
526	USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
527	!!
528	INTEGER, INTENT(in) :: kt ! ocean time-step index
529	!!
530	INTEGER :: ji, jj, jk ! dummy loop indices
531	REAL(wp) :: zcoef0, zuap, zvap ! temporary scalars
532	REAL(wp) :: zalph , zbeta ! " "
533	!!----------------------------------------------------------------------
534
535	IF( kt == nit000 ) THEN
536	IF(lwp) WRITE(numout,*)
537	IF(lwp) WRITE(numout,*) 'dyn:hpg_wdj : hydrostatic pressure gradient trend'
538	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ Weighted Density Jacobian'
539	ENDIF
540
541	! Local constant initialization
542	zcoef0 = - grav * 0.5_wp
543	zalph = 0.5_wp - rn_gamma ! weighting coefficients (alpha=0.5-rn_gamma
544	zbeta = 0.5_wp + rn_gamma ! (beta =1-alpha=0.5+rn_gamma
545
546	! Surface value (no ponderation)
547	DO jj = 2, jpjm1
548	DO ji = fs_2, fs_jpim1 ! vector opt.
549	! hydrostatic pressure gradient along s-surfaces
550	zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * rhd(ji+1,jj ,1) &
551	& - fse3w(ji ,jj ,1) * rhd(ji ,jj ,1) )
552	zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * rhd(ji ,jj+1,1) &
553	& - fse3w(ji ,jj ,1) * rhd(ji, jj ,1) )
554	! s-coordinate pressure gradient correction
555	zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) &
556	& * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
557	zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) &
558	& * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
559	! add to the general momentum trend
560	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
561	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
562	END DO
563	END DO
564
565	! Interior value (2=<jk=<jpkm1) (weighted with zalph & zbeta)
566	DO jk = 2, jpkm1
567	DO jj = 2, jpjm1
568	DO ji = fs_2, fs_jpim1 ! vector opt.
569	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) &
570	& * ( ( fsde3w(ji+1,jj,jk ) + fsde3w(ji,jj,jk ) &
571	& - fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) ) &
572	& * ( zalph * ( rhd (ji+1,jj,jk-1) - rhd (ji,jj,jk-1) ) &
573	& + zbeta * ( rhd (ji+1,jj,jk ) - rhd (ji,jj,jk ) ) ) &
574	& - ( rhd (ji+1,jj,jk ) + rhd (ji,jj,jk ) &
575	& - rhd (ji+1,jj,jk-1) - rhd (ji,jj,jk-1) ) &
576	& * ( zalph * ( fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) ) &
577	& + zbeta * ( fsde3w(ji+1,jj,jk ) - fsde3w(ji,jj,jk ) ) ) )
578	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) &
579	& * ( ( fsde3w(ji,jj+1,jk ) + fsde3w(ji,jj,jk ) &
580	& - fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) ) &
581	& * ( zalph * ( rhd (ji,jj+1,jk-1) - rhd (ji,jj,jk-1) ) &
582	& + zbeta * ( rhd (ji,jj+1,jk ) - rhd (ji,jj,jk ) ) ) &
583	& - ( rhd (ji,jj+1,jk ) + rhd (ji,jj,jk ) &
584	& - rhd (ji,jj+1,jk-1) - rhd (ji,jj,jk-1) ) &
585	& * ( zalph * ( fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) ) &
586	& + zbeta * ( fsde3w(ji,jj+1,jk ) - fsde3w(ji,jj,jk ) ) ) )
587	! add to the general momentum trend
588	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
589	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
590	END DO
591	END DO
592	END DO
593	!
594	END SUBROUTINE hpg_wdj
595
596
597	SUBROUTINE hpg_djc( kt )
598	!!---------------------------------------------------------------------
599	!! * ROUTINE hpg_djc *
600	!!
601	!! ** Method : Density Jacobian with Cubic polynomial scheme
602	!!
603	!! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003
604	!!----------------------------------------------------------------------
605	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
606	USE oce , ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
607	USE wrk_nemo, ONLY: drhox => wrk_3d_1 , dzx => wrk_3d_2
608	USE wrk_nemo, ONLY: drhou => wrk_3d_3 , dzu => wrk_3d_4 , rho_i => wrk_3d_5
609	USE wrk_nemo, ONLY: drhoy => wrk_3d_6 , dzy => wrk_3d_7
610	USE wrk_nemo, ONLY: drhov => wrk_3d_8 , dzv => wrk_3d_9 , rho_j => wrk_3d_10
611	USE wrk_nemo, ONLY: drhoz => wrk_3d_11 , dzz => wrk_3d_12
612	USE wrk_nemo, ONLY: drhow => wrk_3d_13 , dzw => wrk_3d_14
613	USE wrk_nemo, ONLY: rho_k => wrk_3d_15
614	!!
615	INTEGER, INTENT(in) :: kt ! ocean time-step index
616	!!
617	INTEGER :: ji, jj, jk ! dummy loop indices
618	REAL(wp) :: zcoef0, zep, cffw ! temporary scalars
619	REAL(wp) :: z1_10, cffu, cffx ! " "
620	REAL(wp) :: z1_12, cffv, cffy ! " "
621	!!----------------------------------------------------------------------
622
623	IF( wrk_in_use(3, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15) ) THEN
624	CALL ctl_stop('dyn:hpg_djc: requested workspace arrays unavailable') ; RETURN
625	ENDIF
626
627	IF( kt == nit000 ) THEN
628	IF(lwp) WRITE(numout,*)
629	IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend'
630	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme'
631	ENDIF
632
633	! Local constant initialization
634	zcoef0 = - grav * 0.5_wp
635	z1_10 = 1._wp / 10._wp
636	z1_12 = 1._wp / 12._wp
637
638	!----------------------------------------------------------------------------------------
639	! compute and store in provisional arrays elementary vertical and horizontal differences
640	!----------------------------------------------------------------------------------------
641
642	!!bug gm Not a true bug, but... dzz=e3w for dzx, dzy verify what it is really
643
644	DO jk = 2, jpkm1
645	DO jj = 2, jpjm1
646	DO ji = fs_2, fs_jpim1 ! vector opt.
647	drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1)
648	dzz (ji,jj,jk) = fsde3w(ji ,jj ,jk) - fsde3w(ji,jj,jk-1)
649	drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk )
650	dzx (ji,jj,jk) = fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk )
651	drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk )
652	dzy (ji,jj,jk) = fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk )
653	END DO
654	END DO
655	END DO
656
657	!-------------------------------------------------------------------------
658	! compute harmonic averages using eq. 5.18
659	!-------------------------------------------------------------------------
660	zep = 1.e-15
661
662	!!bug gm drhoz not defined at level 1 and used (jk-1 with jk=2)
663	!!bug gm idem for drhox, drhoy et ji=jpi and jj=jpj
664
665	DO jk = 2, jpkm1
666	DO jj = 2, jpjm1
667	DO ji = fs_2, fs_jpim1 ! vector opt.
668	cffw = 2._wp * drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1)
669
670	cffu = 2._wp * drhox(ji+1,jj ,jk) * drhox(ji,jj,jk )
671	cffx = 2._wp * dzx (ji+1,jj ,jk) * dzx (ji,jj,jk )
672
673	cffv = 2._wp * drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk )
674	cffy = 2._wp * dzy (ji ,jj+1,jk) * dzy (ji,jj,jk )
675
676	IF( cffw > zep) THEN
677	drhow(ji,jj,jk) = 2._wp * drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1) &
678	& / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) )
679	ELSE
680	drhow(ji,jj,jk) = 0._wp
681	ENDIF
682
683	dzw(ji,jj,jk) = 2._wp * dzz(ji,jj,jk) * dzz(ji,jj,jk-1) &
684	& / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) )
685
686	IF( cffu > zep ) THEN
687	drhou(ji,jj,jk) = 2._wp * drhox(ji+1,jj,jk) * drhox(ji,jj,jk) &
688	& / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) )
689	ELSE
690	drhou(ji,jj,jk ) = 0._wp
691	ENDIF
692
693	IF( cffx > zep ) THEN
694	dzu(ji,jj,jk) = 2._wp * dzx(ji+1,jj,jk) * dzx(ji,jj,jk) &
695	& / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) )
696	ELSE
697	dzu(ji,jj,jk) = 0._wp
698	ENDIF
699
700	IF( cffv > zep ) THEN
701	drhov(ji,jj,jk) = 2._wp * drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk) &
702	& / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) )
703	ELSE
704	drhov(ji,jj,jk) = 0._wp
705	ENDIF
706
707	IF( cffy > zep ) THEN
708	dzv(ji,jj,jk) = 2._wp * dzy(ji,jj+1,jk) * dzy(ji,jj,jk) &
709	& / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) )
710	ELSE
711	dzv(ji,jj,jk) = 0._wp
712	ENDIF
713
714	END DO
715	END DO
716	END DO
717
718	!----------------------------------------------------------------------------------
719	! apply boundary conditions at top and bottom using 5.36-5.37
720	!----------------------------------------------------------------------------------
721	drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:, 1 ) ) - 0.5_wp * drhow(:,:, 2 )
722	drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:, 1 ) ) - 0.5_wp * drhou(:,:, 2 )
723	drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:, 1 ) ) - 0.5_wp * drhov(:,:, 2 )
724
725	drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1)
726	drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1)
727	drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1)
728
729
730	!--------------------------------------------------------------
731	! Upper half of top-most grid box, compute and store
732	!-------------------------------------------------------------
733
734	!!bug gm : e3w-de3w = 0.5*e3w .... and de3w(2)-de3w(1)=e3w(2) .... to be verified
735	! true if de3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be
736
737	DO jj = 2, jpjm1
738	DO ji = fs_2, fs_jpim1 ! vector opt.
739	rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) ) &
740	& * ( rhd(ji,jj,1) &
741	& + 0.5_wp * ( rhd(ji,jj,2) - rhd(ji,jj,1) ) &
742	& * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) ) &
743	& / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) ) )
744	END DO
745	END DO
746
747	!!bug gm : here also, simplification is possible
748	!!bug gm : optimisation: 1/10 and 1/12 the division should be done before the loop
749
750	DO jk = 2, jpkm1
751	DO jj = 2, jpjm1
752	DO ji = fs_2, fs_jpim1 ! vector opt.
753
754	rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) &
755	& * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) ) &
756	& - grav * z1_10 * ( &
757	& ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) &
758	& * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) &
759	& - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) &
760	& * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) &
761	& )
762
763	rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) &
764	& * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) ) &
765	& - grav* z1_10 * ( &
766	& ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) &
767	& * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) &
768	& - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) &
769	& * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) &
770	& )
771
772	rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) &
773	& * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) ) &
774	& - grav* z1_10 * ( &
775	& ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) &
776	& * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) &
777	& - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) &
778	& * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) &
779	& )
780
781	END DO
782	END DO
783	END DO
784	CALL lbc_lnk(rho_k,'W',1.)
785	CALL lbc_lnk(rho_i,'U',1.)
786	CALL lbc_lnk(rho_j,'V',1.)
787
788
789	! ---------------
790	! Surface value
791	! ---------------
792	DO jj = 2, jpjm1
793	DO ji = fs_2, fs_jpim1 ! vector opt.
794	zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj)
795	zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj)
796	! add to the general momentum trend
797	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
798	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
799	END DO
800	END DO
801
802	! ----------------
803	! interior value (2=<jk=<jpkm1)
804	! ----------------
805	DO jk = 2, jpkm1
806	DO jj = 2, jpjm1
807	DO ji = fs_2, fs_jpim1 ! vector opt.
808	! hydrostatic pressure gradient along s-surfaces
809	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &
810	& + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) ) &
811	& - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) / e1u(ji,jj)
812	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &
813	& + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) &
814	& -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) / e2v(ji,jj)
815	! add to the general momentum trend
816	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
817	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
818	END DO
819	END DO
820	END DO
821	!
822	IF( wrk_not_released(3, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15) ) &
823	CALL ctl_stop('dyn:hpg_djc: failed to release workspace arrays')
824	!
825	END SUBROUTINE hpg_djc
826
827
828	SUBROUTINE hpg_rot( kt )
829	!!---------------------------------------------------------------------
830	!! * ROUTINE hpg_rot *
831	!!
832	!! ** Method : rotated axes scheme (Thiem and Berntsen 2005)
833	!!
834	!! Reference: Thiem & Berntsen, Ocean Modelling, In press, 2005.
835	!!----------------------------------------------------------------------
836	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
837	USE oce , ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
838	USE wrk_nemo, ONLY: zdistr => wrk_2d_1 , zsina => wrk_2d_2 , zcosa => wrk_2d_3
839	USE wrk_nemo, ONLY: zhpiorg => wrk_3d_1 , zhpirot => wrk_3d_2
840	USE wrk_nemo, ONLY: zhpitra => wrk_3d_3 , zhpine => wrk_3d_4
841	USE wrk_nemo, ONLY: zhpjorg => wrk_3d_5 , zhpjrot => wrk_3d_6
842	USE wrk_nemo, ONLY: zhpjtra => wrk_3d_7 , zhpjne => wrk_3d_8
843	!!
844	INTEGER, INTENT(in) :: kt ! ocean time-step index
845	!!
846	INTEGER :: ji, jj, jk ! dummy loop indices
847	REAL(wp) :: zforg, zcoef0, zuap, zmskd1, zmskd1m ! temporary scalar
848	REAL(wp) :: zfrot , zvap, zmskd2, zmskd2m ! " "
849	!!----------------------------------------------------------------------
850
851	IF( wrk_in_use(2, 1,2,3) .OR. &
852	wrk_in_use(3, 1,2,3,4,5,6,7,8) ) THEN
853	CALL ctl_stop('dyn:hpg_rot: requested workspace arrays unavailable') ; RETURN
854	ENDIF
855
856	IF( kt == nit000 ) THEN
857	IF(lwp) WRITE(numout,*)
858	IF(lwp) WRITE(numout,*) 'dyn:hpg_rot : hydrostatic pressure gradient trend'
859	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, rotated axes scheme used'
860	ENDIF
861
862	! -------------------------------
863	! Local constant initialization
864	! -------------------------------
865	zcoef0 = - grav * 0.5_wp
866	zforg = 0.95_wp
867	zfrot = 1._wp - zforg
868
869	! inverse of the distance between 2 diagonal T-points (defined at F-point) (here zcoef0/distance)
870	zdistr(:,:) = zcoef0 / SQRT( e1f(:,:)e1f(:,:) + e2f(:,:)e1f(:,:) )
871
872	! sinus and cosinus of diagonal angle at F-point
873	zsina(:,:) = ATAN2( e2f(:,:), e1f(:,:) )
874	zcosa(:,:) = COS( zsina(:,:) )
875	zsina(:,:) = SIN( zsina(:,:) )
876
877	! ---------------
878	! Surface value
879	! ---------------
880	! compute and add to the general trend the pressure gradients along the axes
881	DO jj = 2, jpjm1
882	DO ji = fs_2, fs_jpim1 ! vector opt.
883	! hydrostatic pressure gradient along s-surfaces
884	zhpiorg(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,1) * rhd(ji+1,jj,1) &
885	& - fse3t(ji ,jj,1) * rhd(ji ,jj,1) )
886	zhpjorg(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,1) * rhd(ji,jj+1,1) &
887	& - fse3t(ji,jj ,1) * rhd(ji,jj ,1) )
888	! s-coordinate pressure gradient correction
889	zuap = -zcoef0 * ( rhd (ji+1,jj ,1) + rhd (ji,jj,1) ) &
890	& * ( fsdept(ji+1,jj ,1) - fsdept(ji,jj,1) ) / e1u(ji,jj)
891	zvap = -zcoef0 * ( rhd (ji ,jj+1,1) + rhd (ji,jj,1) ) &
892	& * ( fsdept(ji ,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj)
893	! add to the general momentum trend
894	ua(ji,jj,1) = ua(ji,jj,1) + zforg * ( zhpiorg(ji,jj,1) + zuap )
895	va(ji,jj,1) = va(ji,jj,1) + zforg * ( zhpjorg(ji,jj,1) + zvap )
896	END DO
897	END DO
898
899	! compute the pressure gradients in the diagonal directions
900	DO jj = 1, jpjm1
901	DO ji = 1, fs_jpim1 ! vector opt.
902	zmskd1 = tmask(ji+1,jj+1,1) * tmask(ji ,jj,1) ! mask in the 1st diagnonal
903	zmskd2 = tmask(ji ,jj+1,1) * tmask(ji+1,jj,1) ! mask in the 2nd diagnonal
904	! hydrostatic pressure gradient along s-surfaces
905	zhpitra(ji,jj,1) = zdistr(ji,jj) * zmskd1 * ( fse3t(ji+1,jj+1,1) * rhd(ji+1,jj+1,1) &
906	& - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) )
907	zhpjtra(ji,jj,1) = zdistr(ji,jj) * zmskd2 * ( fse3t(ji ,jj+1,1) * rhd(ji ,jj+1,1) &
908	& - fse3t(ji+1,jj ,1) * rhd(ji+1,jj ,1) )
909	! s-coordinate pressure gradient correction
910	zuap = -zdistr(ji,jj) * zmskd1 * ( rhd (ji+1,jj+1,1) + rhd (ji ,jj,1) ) &
911	& * ( fsdept(ji+1,jj+1,1) - fsdept(ji ,jj,1) )
912	zvap = -zdistr(ji,jj) * zmskd2 * ( rhd (ji ,jj+1,1) + rhd (ji+1,jj,1) ) &
913	& * ( fsdept(ji ,jj+1,1) - fsdept(ji+1,jj,1) )
914	! back rotation
915	zhpine(ji,jj,1) = zcosa(ji,jj) * ( zhpitra(ji,jj,1) + zuap ) &
916	& - zsina(ji,jj) * ( zhpjtra(ji,jj,1) + zvap )
917	zhpjne(ji,jj,1) = zsina(ji,jj) * ( zhpitra(ji,jj,1) + zuap ) &
918	& + zcosa(ji,jj) * ( zhpjtra(ji,jj,1) + zvap )
919	END DO
920	END DO
921
922	! interpolate and add to the general trend the diagonal gradient
923	DO jj = 2, jpjm1
924	DO ji = fs_2, fs_jpim1 ! vector opt.
925	! averaging
926	zhpirot(ji,jj,1) = 0.5 * ( zhpine(ji,jj,1) + zhpine(ji ,jj-1,1) )
927	zhpjrot(ji,jj,1) = 0.5 * ( zhpjne(ji,jj,1) + zhpjne(ji-1,jj ,1) )
928	! add to the general momentum trend
929	ua(ji,jj,1) = ua(ji,jj,1) + zfrot * zhpirot(ji,jj,1)
930	va(ji,jj,1) = va(ji,jj,1) + zfrot * zhpjrot(ji,jj,1)
931	END DO
932	END DO
933
934	! -----------------
935	! 2. interior value (2=<jk=<jpkm1)
936	! -----------------
937	! compute and add to the general trend the pressure gradients along the axes
938	DO jk = 2, jpkm1
939	DO jj = 2, jpjm1
940	DO ji = fs_2, fs_jpim1 ! vector opt.
941	! hydrostatic pressure gradient along s-surfaces
942	zhpiorg(ji,jj,jk) = zhpiorg(ji,jj,jk-1) &
943	& + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk ) * rhd(ji+1,jj,jk ) &
944	& - fse3t(ji ,jj,jk ) * rhd(ji ,jj,jk ) &
945	& + fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1) &
946	& - fse3t(ji ,jj,jk-1) * rhd(ji ,jj,jk-1) )
947	zhpjorg(ji,jj,jk) = zhpjorg(ji,jj,jk-1) &
948	& + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk ) &
949	& - fse3t(ji,jj ,jk ) * rhd(ji,jj, jk ) &
950	& + fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) &
951	& - fse3t(ji,jj ,jk-1) * rhd(ji,jj, jk-1) )
952	! s-coordinate pressure gradient correction
953	zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) ) &
954	& * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)
955	zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) ) &
956	& * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)
957	! add to the general momentum trend
958	ua(ji,jj,jk) = ua(ji,jj,jk) + zforg*( zhpiorg(ji,jj,jk) + zuap )
959	va(ji,jj,jk) = va(ji,jj,jk) + zforg*( zhpjorg(ji,jj,jk) + zvap )
960	END DO
961	END DO
962	END DO
963
964	! compute the pressure gradients in the diagonal directions
965	DO jk = 2, jpkm1
966	DO jj = 1, jpjm1
967	DO ji = 1, fs_jpim1 ! vector opt.
968	zmskd1 = tmask(ji+1,jj+1,jk ) * tmask(ji ,jj,jk ) ! level jk mask in the 1st diagnonal
969	zmskd1m = tmask(ji+1,jj+1,jk-1) * tmask(ji ,jj,jk-1) ! level jk-1 " "
970	zmskd2 = tmask(ji ,jj+1,jk ) * tmask(ji+1,jj,jk ) ! level jk mask in the 2nd diagnonal
971	zmskd2m = tmask(ji ,jj+1,jk-1) * tmask(ji+1,jj,jk-1) ! level jk-1 " "
972	! hydrostatic pressure gradient along s-surfaces
973	zhpitra(ji,jj,jk) = zhpitra(ji,jj,jk-1) &
974	& + zdistr(ji,jj) * zmskd1 * ( fse3t(ji+1,jj+1,jk ) * rhd(ji+1,jj+1,jk) &
975	& -fse3t(ji ,jj ,jk ) * rhd(ji ,jj ,jk) ) &
976	& + zdistr(ji,jj) * zmskd1m * ( fse3t(ji+1,jj+1,jk-1) * rhd(ji+1,jj+1,jk-1) &
977	& -fse3t(ji ,jj ,jk-1) * rhd(ji ,jj ,jk-1) )
978	zhpjtra(ji,jj,jk) = zhpjtra(ji,jj,jk-1) &
979	& + zdistr(ji,jj) * zmskd2 * ( fse3t(ji ,jj+1,jk ) * rhd(ji ,jj+1,jk) &
980	& -fse3t(ji+1,jj ,jk ) * rhd(ji+1,jj, jk) ) &
981	& + zdistr(ji,jj) * zmskd2m * ( fse3t(ji ,jj+1,jk-1) * rhd(ji ,jj+1,jk-1) &
982	& -fse3t(ji+1,jj ,jk-1) * rhd(ji+1,jj, jk-1) )
983	! s-coordinate pressure gradient correction
984	zuap = - zdistr(ji,jj) * zmskd1 * ( rhd (ji+1,jj+1,jk) + rhd (ji ,jj,jk) ) &
985	& * ( fsdept(ji+1,jj+1,jk) - fsdept(ji ,jj,jk) )
986	zvap = - zdistr(ji,jj) * zmskd2 * ( rhd (ji ,jj+1,jk) + rhd (ji+1,jj,jk) ) &
987	& * ( fsdept(ji ,jj+1,jk) - fsdept(ji+1,jj,jk) )
988	! back rotation
989	zhpine(ji,jj,jk) = zcosa(ji,jj) * ( zhpitra(ji,jj,jk) + zuap ) &
990	& - zsina(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap )
991	zhpjne(ji,jj,jk) = zsina(ji,jj) * ( zhpitra(ji,jj,jk) + zuap ) &
992	& + zcosa(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap )
993	END DO
994	END DO
995	END DO
996
997	! interpolate and add to the general trend
998	DO jk = 2, jpkm1
999	DO jj = 2, jpjm1
1000	DO ji = fs_2, fs_jpim1 ! vector opt.
1001	! averaging
1002	zhpirot(ji,jj,jk) = 0.5 * ( zhpine(ji,jj,jk) + zhpine(ji ,jj-1,jk) )
1003	zhpjrot(ji,jj,jk) = 0.5 * ( zhpjne(ji,jj,jk) + zhpjne(ji-1,jj ,jk) )
1004	! add to the general momentum trend
1005	ua(ji,jj,jk) = ua(ji,jj,jk) + zfrot * zhpirot(ji,jj,jk)
1006	va(ji,jj,jk) = va(ji,jj,jk) + zfrot * zhpjrot(ji,jj,jk)
1007	END DO
1008	END DO
1009	END DO
1010	!
1011	IF( wrk_not_released(2, 1,2,3) .OR. &
1012	wrk_not_released(3, 1,2,3,4,5,6,7,8) ) CALL ctl_stop('dyn:hpg_rot: failed to release workspace arrays')
1013	!
1014	END SUBROUTINE hpg_rot
1015
1016
1017	SUBROUTINE hpg_prj( kt )
1018	!!---------------------------------------------------------------------
1019	!! * ROUTINE hpg_prj *
1020	!!
1021	!! ** Method : s-coordinate case.
1022	!! Reformulate the horizontal hydrostatical pressure gradient
1023	!! term using Pressure Jacobian. A new correction term
1024	!! is developed to eliminate the sigma-coordinate error.
1025	!!
1026	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
1027	!! - Save the trend (l_trddyn=T)
1028	!!
1029	!!----------------------------------------------------------------------
1030
1031	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
1032	USE oce , ONLY: fsde3w_tmp => ta ! (ta,sa) used as 3D workspace
1033	USE oce , ONLY: rhd_tmp => sa
1034	USE wrk_nemo, ONLY: zhpi => wrk_3d_3
1035	USE wrk_nemo, ONLY: zu => wrk_3d_4
1036	USE wrk_nemo, ONLY: zv => wrk_3d_5
1037	USE wrk_nemo, ONLY: fsp => wrk_3d_6
1038	USE wrk_nemo, ONLY: xsp => wrk_3d_7
1039	USE wrk_nemo, ONLY: asp => wrk_3d_8
1040	USE wrk_nemo, ONLY: bsp => wrk_3d_9
1041	USE wrk_nemo, ONLY: csp => wrk_3d_10
1042	USE wrk_nemo, ONLY: dsp => wrk_3d_11
1043	!!
1044	!!----------------------------------------------------------------------
1045	!!
1046	INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear
1047	INTEGER, INTENT(in) :: kt ! ocean time-step index
1048	!!
1049	INTEGER :: ji, jj, jk, jkk ! dummy loop indices
1050	REAL(wp) :: zcoef0, znad ! temporary scalars
1051	!!
1052	!! The local varialbes for the correction term
1053	INTEGER :: jke, jkw, jkn, jks, jk1, jis, jid, jjs, jjd
1054	REAL(wp) :: zze, zzw, zzn, zzs, rre, rrw, rrn, rrs
1055	REAL(wp) :: zuijk, zvijk, pe, pw, pwes, pwed, pn, ps, pnss, pnsd, deps
1056	REAL(wp) :: rhde, rhdw, rhdn, rhds,rhdt1, rhdt2
1057	REAL(wp) :: dpdx1, dpdx2, dpdy1, dpdy2
1058	INTEGER :: bhitwe, bhitns
1059	!!----------------------------------------------------------------------
1060
1061	IF( wrk_in_use(3, 3,4,5,6,7,8,9,10,11) ) THEN
1062	CALL ctl_stop('dyn:hpg_prj: requested workspace arrays unavailable') ; RETURN
1063	ENDIF
1064
1065	IF( kt == nit000 ) THEN
1066	IF(lwp) WRITE(numout,*)
1067	IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend'
1068	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, cubic spline pressure Jacobian'
1069	ENDIF
1070
1071	!!----------------------------------------------------------------------
1072	! Local constant initialization
1073	zcoef0 = - grav
1074	znad = 0.0_wp
1075	IF(lk_vvl) znad = 1._wp
1076
1077	! Save fsde3w and rhd
1078	fsde3w_tmp(:,:,:) = fsde3w(:,:,:)
1079	rhd_tmp(:,:,:) = rhd(:,:,:)
1080
1081	! Clean 3-D work arraies
1082	zhpi(:,:,:) = 0.
1083
1084	!how to add vector opt.? N.B., jpi&jpi rather than jpim1&jpjm1 are needed here
1085
1086	! Preparing vertical density profile for hybrid-sco coordinate
1087	DO jj = 1, jpj
1088	DO ji = 1, jpi
1089	jk = mbathy(ji,jj)
1090	IF(jk <= 0) THEN; rhd(ji,jj,:) = 0._wp
1091	ELSE IF(jk == 1) THEN; rhd(ji,jj, jk+1:jpk) = rhd(ji,jj,jk)
1092	ELSE IF(jk < jpkm1) THEN
1093	DO jkk = jk+1, jpk
1094	rhd(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1),&
1095	fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2))
1096	END DO
1097	END IF
1098	END DO
1099	END DO
1100
1101	DO jj = 1, jpj
1102	DO ji = 1, jpi
1103	fsde3w(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1)
1104	fsde3w(ji,jj,1) = fsde3w(ji,jj,1) - sshn(ji,jj) * znad
1105	DO jk = 2, jpk
1106	fsde3w(ji,jj,jk) = fsde3w(ji,jj,jk-1) + fse3w(ji,jj,jk)
1107	END DO
1108	END DO
1109	END DO
1110
1111	DO jk = 1, jpkm1
1112	DO jj = 1, jpj
1113	DO ji = 1, jpi
1114	fsp(ji,jj,jk) = rhd(ji,jj,jk)
1115	xsp(ji,jj,jk) = fsde3w(ji,jj,jk)
1116	END DO
1117	END DO
1118	END DO
1119
1120	! ! Construct the vertical density profile with the
1121	! !Constrained cubic spline interpolation
1122	CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type)
1123
1124	! Calculate the hydrostatic pressure at T(ji,jj,1)
1125	DO jj = 2, jpj
1126	DO ji = 2, jpi
1127	rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), &
1128	bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) &
1129	* 0.5_wp * fsde3w(ji,jj,1)
1130	rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water
1131
1132	! ! assuming linear profile across the top half surface layer
1133	zhpi(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) * rhdt1
1134	END DO
1135	END DO
1136
1137	! Calculate the pressure at T(ji,jj,2:jpkm1)
1138	DO jk = 2, jpkm1
1139	DO jj = 2, jpj
1140	DO ji = 2, jpi
1141	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + &
1142	integ2(fsde3w(ji,jj,jk-1),fsde3w(ji,jj,jk),&
1143	asp(ji,jj,jk-1),bsp(ji,jj,jk-1),&
1144	csp(ji,jj,jk-1),dsp(ji,jj,jk-1))
1145	END DO
1146	END DO
1147	END DO
1148
1149	! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1)
1150
1151	DO jj = 2, jpjm1
1152	DO ji = 2, jpim1
1153	zu(ji,jj,1) = - ( 0.5_wp * fse3uw(ji,jj,1) - sshu_n(ji,jj) * znad)
1154	zv(ji,jj,1) = - ( 0.5_wp * fse3vw(ji,jj,1) - sshv_n(ji,jj) * znad)
1155	END DO
1156	END DO
1157
1158	DO jk = 2, jpkm1
1159	DO jj = 2, jpjm1
1160	DO ji = 2, jpim1
1161	zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3uw(ji,jj,jk)
1162	zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3vw(ji,jj,jk)
1163	END DO
1164	END DO
1165	END DO
1166
1167	! Start pressure integration
1168
1169	DO jk = 1, jpkm1
1170	DO jj = 2, jpjm1
1171	DO ji = 2, jpim1
1172	pwes = 0._wp; pwed = 0._wp
1173	pnss = 0._wp; pnsd = 0._wp
1174	zuijk = zu(ji,jj,jk)
1175	zvijk = zv(ji,jj,jk)
1176
1177	!!!!! for u equation
1178	IF( -fsde3w(ji+1,jj,mbathy(ji+1,jj)) >= -fsde3w(ji,jj,mbathy(ji,jj)) ) THEN
1179	jis = ji + 1; jid = ji
1180	ELSE
1181	jis = ji; jid = ji +1
1182	ENDIF
1183
1184	! !integrate the pressure on the shallow side
1185	jk1 = jk
1186	bhitwe = 0
1187	IF( jk1 == mbathy(jis,jj) ) THEN
1188	bhitwe = 1
1189	ELSE
1190	DO WHILE ( -fsde3w(jis,jj,jk1+1) > zuijk )
1191	pwes = pwes + &
1192	integ2(fsde3w(jis,jj,jk1), fsde3w(jis,jj,jk1+1),&
1193	asp(jis,jj,jk1) , bsp(jis,jj,jk1), &
1194	csp(jis,jj,jk1) , dsp(jis,jj,jk1))
1195	jk1 = jk1 + 1
1196	IF( jk1 == mbathy(jis,jj) ) THEN
1197	bhitwe = 1
1198	EXIT
1199	END IF
1200	END DO
1201	ENDIF
1202
1203	IF( bhitwe /= 1 ) THEN
1204	pwes = pwes + &
1205	integ2(fsde3w(jis,jj,jk1), -zuijk, &
1206	asp(jis,jj,jk1), bsp(jis,jj,jk1),&
1207	csp(jis,jj,jk1), dsp(jis,jj,jk1))
1208	ELSE
1209	zuijk = -fsde3w(jis,jj,jk1)
1210	ENDIF
1211
1212	! !integrate the pressure on the deep side
1213	jk1 = jk
1214	bhitwe = 0
1215	IF( jk1 == 1 ) THEN
1216	bhitwe = 1
1217	ELSE
1218	DO WHILE ( -fsde3w(jid,jj,jk1-1) < zuijk )
1219	pwed = pwed + &
1220	integ2(fsde3w(jid,jj,jk1-1), fsde3w(jid,jj,jk1),&
1221	asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), &
1222	csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1))
1223	jk1 = jk1 - 1
1224	IF( jk1 == 1 ) THEN
1225	bhitwe = 1
1226	EXIT
1227	END IF
1228	END DO
1229	ENDIF
1230
1231	IF( bhitwe /= 1 ) THEN
1232	pwed = pwed + &
1233	integ2(-zuijk, fsde3w(jid,jj,jk1),&
1234	asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), &
1235	csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1))
1236	ELSE
1237	deps = fsde3w(jid,jj,1) + min(zuijk, sshn(jid,jj)*znad)
1238	rhdt1 = rhd(jid,jj,1) - interp3(fsde3w(jid,jj,1), asp(jid,jj,1), &
1239	bsp(jid,jj,1), csp(jid,jj,1), &
1240	dsp(jid,jj,1) ) * deps
1241	rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water
1242	pwed = pwed + 0.5_wp * (rhd(jid,jj,1) + rhdt1) * deps
1243	ENDIF
1244
1245	IF( jid > jis ) THEN
1246	pe = pwed; pw = pwes
1247	ELSE
1248	pe = pwes; pw = pwed
1249	ENDIF
1250
1251	dpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk))
1252	IF( lk_vvl ) THEN
1253	dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe + (sshn(ji+1,jj)-sshn(ji,jj)))
1254	ELSE
1255	dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe)
1256	ENDIF
1257
1258	ua(ji,jj,jk) = ua(ji,jj,jk) + (dpdx1 + dpdx2) * &
1259	& umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk)
1260
1261	!!!!! for v equation
1262
1263	IF( -fsde3w(ji,jj+1,mbathy(ji,jj+1)) >= -fsde3w(ji,jj,mbathy(ji,jj)) ) THEN
1264	jjs = jj + 1; jjd = jj
1265	ELSE
1266	jjs = jj ; jjd = jj + 1
1267	ENDIF
1268
1269	! !integrate the pressure on the shallow side
1270	jk1 = jk
1271	bhitns = 0
1272	IF( jk1 == mbathy(ji,jjs) ) THEN
1273	bhitns = 1
1274	ELSE
1275	DO WHILE ( -fsde3w(ji,jjs,jk1+1) > zvijk )
1276	pnss = pnss + &
1277	integ2(fsde3w(ji,jjs,jk1), fsde3w(ji,jjs,jk1+1),&
1278	asp(ji,jjs,jk1), bsp(ji,jjs,jk1), &
1279	csp(ji,jjs,jk1), dsp(ji,jjs,jk1) )
1280	jk1 = jk1 + 1
1281	IF( jk1 == mbathy(ji,jjs) ) THEN
1282	bhitns = 1
1283	EXIT
1284	END IF
1285	END DO
1286	ENDIF
1287
1288	IF( bhitns /= 1 ) THEN
1289	pnss = pnss + &
1290	integ2(fsde3w(ji,jjs,jk1), -zvijk, &
1291	asp(ji,jjs,jk1), bsp(ji,jjs,jk1), &
1292	csp(ji,jjs,jk1), dsp(ji,jjs,jk1) )
1293	ELSE
1294	zvijk = -fsde3w(ji,jjs,jk1)
1295	ENDIF
1296
1297	! !integrate the pressure on the deep side
1298	jk1 = jk
1299	bhitns = 0
1300	IF( jk1 == 1 ) THEN
1301	bhitns = 1
1302	ELSE
1303	DO WHILE ( -fsde3w(ji,jjd,jk1-1) < zvijk )
1304	pnsd = pnsd + &
1305	integ2(fsde3w(ji,jjd,jk1-1), fsde3w(ji,jjd,jk1), &
1306	asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), &
1307	csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
1308	jk1 = jk1 - 1
1309	IF( jk1 == 1 ) THEN
1310	bhitns = 1
1311	EXIT
1312	END IF
1313	END DO
1314	ENDIF
1315
1316	IF( bhitns /= 1 ) THEN
1317	pnsd = pnsd + &
1318	integ2(-zvijk, fsde3w(ji,jjd,jk1), &
1319	asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), &
1320	csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )
1321	ELSE
1322	deps = fsde3w(ji,jjd,1) + min(zvijk, sshn(ji,jjd)*znad)
1323	rhdt1 = rhd(ji,jjd,1) - interp3(fsde3w(ji,jjd,1), asp(ji,jjd,1), &
1324	bsp(ji,jjd,1), csp(ji,jjd,1), &
1325	dsp(ji,jjd,1) ) * deps
1326	rhdt1 = max(rhdt1, 1000._wp - rau0) ! no lighter than fresh water
1327	pnsd = pnsd + 0.5_wp * (rhd(ji,jjd,1) + rhdt1) * deps
1328	ENDIF
1329
1330	IF( jjd > jjs ) THEN
1331	pn = pnsd; ps = pnss
1332	ELSE
1333	pn = pnss; ps = pnsd
1334	ENDIF
1335
1336	dpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk))
1337	if( lk_vvl ) then
1338	dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn + (sshn(ji,jj+1)-sshn(ji,jj)))
1339	else
1340	dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn )
1341	end if
1342
1343	va(ji,jj,jk) = va(ji,jj,jk) + (dpdy1 + dpdy2)* &
1344	& vmask(ji,jj,jk)tmask(ji,jj,jk)tmask(ji,jj+1,jk)
1345
1346
1347	END DO
1348	END DO
1349	END DO
1350
1351	! Restore fsde3w and rhd
1352	fsde3w(:,:,:) = fsde3w_tmp(:,:,:)
1353	rhd(:,:,:) = rhd_tmp(:,:,:)
1354
1355	!
1356	IF( wrk_not_released(3, 3,4,5,6,7,8,9,10,11) ) &
1357	CALL ctl_stop('dyn:hpg_prj: failed to release workspace arrays')
1358	!
1359	END SUBROUTINE hpg_prj
1360
1361	SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type)
1362	!!----------------------------------------------------------------------
1363	!! * ROUTINE cspline *
1364	!!
1365	!! ** Purpose : constrained cubic spline interpolation
1366	!!
1367	!! ** Method : f(x) = asp + bspx + cspx^2 + dsp*x^3
1368	!! Reference: K.W. Brodlie, A review of mehtods for curve and function
1369	!! drawing, 1980
1370	!!
1371	!!----------------------------------------------------------------------
1372	IMPLICIT NONE
1373	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: fsp, xsp ! value and coordinate
1374	REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of
1375	! the interpoated function
1376	INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline
1377	! 2: Linear
1378
1379	! Local Variables
1380	INTEGER :: ji, jj, jk ! dummy loop indices
1381	INTEGER :: jpi, jpj, jpkm1
1382	REAL(wp) :: zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp
1383	REAL(wp) :: zdxtmp1, zdxtmp2, zalpha
1384	REAL(wp) :: zdf(size(fsp,3))
1385	!!----------------------------------------------------------------------
1386
1387	jpi = size(fsp,1)
1388	jpj = size(fsp,2)
1389	jpkm1 = size(fsp,3) - 1
1390
1391	! Clean output arrays
1392	asp = 0.0_wp
1393	bsp = 0.0_wp
1394	csp = 0.0_wp
1395	dsp = 0.0_wp
1396
1397
1398	Do ji = 1, jpi
1399	Do jj = 1, jpj
1400
1401	If (polynomial_type == 1) Then !Constrained Cubic Spline
1402	Do jk = 2, jpkm1-1
1403	zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1)
1404	zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1405	zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1
1406	zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2
1407
1408	zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp
1409
1410	If(zdf1 * zdf2 <= 0._wp) Then
1411	zdf(jk) = 0._wp
1412	Else
1413	zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 )
1414	End If
1415	End Do
1416
1417	zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - &
1418	& 0.5_wp * zdf(2)
1419	zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / &
1420	& ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - &
1421	& 0.5_wp * zdf(jpkm1 - 1)
1422
1423	Do jk = 1, jpkm1 - 1
1424	zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1425	ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp
1426	ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp
1427	zddf1 = -2._wp * ztmp1 + ztmp2
1428	ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp
1429	zddf2 = 2._wp * ztmp1 - ztmp2
1430
1431	dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp
1432	csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp
1433	bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - &
1434	& csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - &
1435	& dsp(ji,jj,jk) * ( xsp(ji,jj,jk+1)**2 + &
1436	& xsp(ji,jj,jk+1) * xsp(ji,jj,jk) + &
1437	& xsp(ji,jj,jk)**2 )
1438	asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) - &
1439	& csp(ji,jj,jk) * xsp(ji,jj,jk)**2 - &
1440	& dsp(ji,jj,jk) * xsp(ji,jj,jk)**3
1441	End Do
1442
1443	Else If (polynomial_type == 2) Then !Linear
1444
1445	Do jk = 1, jpkm1-1
1446	zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk)
1447	ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk)
1448
1449	dsp(ji,jj,jk) = 0._wp
1450	csp(ji,jj,jk) = 0._wp
1451	bsp(ji,jj,jk) = ztmp1 / zdxtmp
1452	asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk)
1453	End Do
1454
1455	Else
1456	CALL ctl_stop( 'invalid polynomial type in cspline' )
1457	End If
1458
1459	End Do
1460	End Do
1461
1462	End Subroutine cspline
1463
1464
1465	FUNCTION interp1(x, xl, xr, fl, fr) RESULT(f)
1466	!!----------------------------------------------------------------------
1467	!! * ROUTINE interp1 *
1468	!!
1469	!! ** Purpose : 1-d linear interpolation
1470	!!
1471	!! ** Method :
1472	!! interpolation is straight forward
1473	!! extrapolation is also permitted (no value limit)
1474	!!
1475	!! H.Liu, Jan 2009, POL
1476	!!----------------------------------------------------------------------
1477	IMPLICIT NONE
1478	REAL(wp), INTENT(in) :: x, xl, xr, fl, fr
1479	REAL(wp) :: f ! result of the interpolation (extrapolation)
1480	REAL(wp) :: zdeltx
1481	!!----------------------------------------------------------------------
1482
1483	zdeltx = xr - xl
1484	IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN
1485	f = 0.5_wp * (fl + fr)
1486	ELSE
1487	f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx
1488	END IF
1489
1490	END FUNCTION interp1
1491
1492	FUNCTION interp2(x, a, b, c, d) RESULT(f)
1493	!!----------------------------------------------------------------------
1494	!! * ROUTINE interp1 *
1495	!!
1496	!! ** Purpose : 1-d constrained cubic spline interpolation
1497	!!
1498	!! ** Method : cubic spline interpolation
1499	!!
1500	!!----------------------------------------------------------------------
1501	IMPLICIT NONE
1502	REAL(wp), INTENT(in) :: x, a, b, c, d
1503	REAL(wp) :: f ! value from the interpolation
1504	!!----------------------------------------------------------------------
1505
1506	f = a + x* ( b + x * ( c + d * x ) )
1507
1508	END FUNCTION interp2
1509
1510
1511	FUNCTION interp3(x, a, b, c, d) RESULT(f)
1512	!!----------------------------------------------------------------------
1513	!! * ROUTINE interp1 *
1514	!!
1515	!! ** Purpose : deriavtive of a cubic spline function
1516	!!
1517	!! ** Method :
1518	!!
1519	!!----------------------------------------------------------------------
1520	IMPLICIT NONE
1521	REAL(wp), INTENT(in) :: x, a, b, c, d
1522	REAL(wp) :: f ! value from the interpolation
1523	!!----------------------------------------------------------------------
1524
1525	f = b + x * ( 2._wp * c + 3._wp * d * x)
1526
1527	END FUNCTION interp3
1528
1529
1530	FUNCTION integ2(xl, xr, a, b, c, d) RESULT(f)
1531	!!----------------------------------------------------------------------
1532	!! * ROUTINE interp1 *
1533	!!
1534	!! ** Purpose : 1-d constrained cubic spline integration
1535	!!
1536	!! ** Method : integrate polynomial a+bx+cx^2+dx^3 from xl to xr
1537	!!
1538	!!----------------------------------------------------------------------
1539	IMPLICIT NONE
1540	REAL(wp), INTENT(in) :: xl, xr, a, b, c, d
1541	REAL(wp) :: za1, za2, za3
1542	REAL(wp) :: f ! integration result
1543	!!----------------------------------------------------------------------
1544
1545	za1 = 0.5_wp * b
1546	za2 = c / 3.0_wp
1547	za3 = 0.25_wp * d
1548
1549	f = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - &
1550	& xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) )
1551
1552	END FUNCTION integ2
1553
1554
1555	!!======================================================================
1556	END MODULE dynhpg

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

Download in other formats: