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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 7698

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

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

Download in other formats: