Context Navigation

dynhpg.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

File size: 71.6 KB

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

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/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 11101

Download in other formats: