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

Context Navigation

source: branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 4747

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

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

Download in other formats: