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

Context Navigation

source: branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 4328

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

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

Download in other formats: