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/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 5212

Last change on this file since 5212 was 5212, checked in by hliu, 9 years ago
A serious mistake in dynhpg: hpg_prj, which should also be corrected in v3.6 main trunk
Property svn:keywords set to `Id`
File size: 59.8 KB

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

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

Download in other formats: