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

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 4400

Last change on this file since 4400 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 56.3 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
	13	!! ! Original code for hpg_ctl, hpg_hel hpg_wdj, hpg_djc, hpg_rot
	14	!! - ! 2005-11 (G. Madec) style & small optimisation
	15	!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
[503]	16	!!----------------------------------------------------------------------
[3]	17
	18	!!----------------------------------------------------------------------
[455]	19	!! dyn_hpg : update the momentum trend with the now horizontal
[3]	20	!! gradient of the hydrostatic pressure
[2528]	21	!! dyn_hpg_init : initialisation and control of options
[455]	22	!! hpg_zco : z-coordinate scheme
	23	!! hpg_zps : z-coordinate plus partial steps (interpolation)
	24	!! hpg_sco : s-coordinate (standard jacobian formulation)
	25	!! hpg_hel : s-coordinate (helsinki modification)
	26	!! hpg_wdj : s-coordinate (weighted density jacobian)
	27	!! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial)
	28	!! hpg_rot : s-coordinate (ROTated axes scheme)
[3]	29	!!----------------------------------------------------------------------
	30	USE oce ! ocean dynamics and tracers
	31	USE dom_oce ! ocean space and time domain
	32	USE phycst ! physical constants
[216]	33	USE trdmod ! ocean dynamics trends
	34	USE trdmod_oce ! ocean variables trends
[2715]	35	USE in_out_manager ! I/O manager
[258]	36	USE prtctl ! Print control
[455]	37	USE lbclnk ! lateral boundary condition
[2715]	38	USE lib_mpp ! MPP library
[3]	39
	40	IMPLICIT NONE
	41	PRIVATE
	42
[2528]	43	PUBLIC dyn_hpg ! routine called by step module
	44	PUBLIC dyn_hpg_init ! routine called by opa module
[3]	45
[1601]	46	! !!* Namelist namdyn_hpg : hydrostatic pressure gradient
	47	LOGICAL , PUBLIC :: ln_hpg_zco = .TRUE. !: z-coordinate - full steps
	48	LOGICAL , PUBLIC :: ln_hpg_zps = .FALSE. !: z-coordinate - partial steps (interpolation)
	49	LOGICAL , PUBLIC :: ln_hpg_sco = .FALSE. !: s-coordinate (standard jacobian formulation)
	50	LOGICAL , PUBLIC :: ln_hpg_hel = .FALSE. !: s-coordinate (helsinki modification)
	51	LOGICAL , PUBLIC :: ln_hpg_wdj = .FALSE. !: s-coordinate (weighted density jacobian)
	52	LOGICAL , PUBLIC :: ln_hpg_djc = .FALSE. !: s-coordinate (Density Jacobian with Cubic polynomial)
	53	LOGICAL , PUBLIC :: ln_hpg_rot = .FALSE. !: s-coordinate (ROTated axes scheme)
[2528]	54	REAL(wp), PUBLIC :: rn_gamma = 0._wp !: weighting coefficient
[1601]	55	LOGICAL , PUBLIC :: ln_dynhpg_imp = .FALSE. !: semi-implicite hpg flag
[455]	56
[1601]	57	INTEGER :: nhpg = 0 ! = 0 to 6, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags)
[455]	58
[3211]	59	!! * Control permutation of array indices
	60	# include "oce_ftrans.h90"
	61	# include "dom_oce_ftrans.h90"
	62
[3]	63	!! * Substitutions
	64	# include "domzgr_substitute.h90"
	65	# include "vectopt_loop_substitute.h90"
	66	!!----------------------------------------------------------------------
[2528]	67	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
	68	!! $Id$
	69	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	70	!!----------------------------------------------------------------------
	71	CONTAINS
	72
	73	SUBROUTINE dyn_hpg( kt )
	74	!!---------------------------------------------------------------------
	75	!! * ROUTINE dyn_hpg *
	76	!!
[455]	77	!! ** Method : Call the hydrostatic pressure gradient routine
[503]	78	!! using the scheme defined in the namelist
[455]	79	!!
	80	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
	81	!! - Save the trend (l_trddyn=T)
[503]	82	!!----------------------------------------------------------------------
[2715]	83	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	84	USE wrk_nemo, ONLY: ztrdu => wrk_3d_1 , ztrdv => wrk_3d_2 ! 3D workspace
[3211]	85	!! DCSE_NEMO: need additional directives for renamed module variables
	86	!FTRANS ztrdu :I :I :z
	87	!FTRANS ztrdv :I :I :z
	88
[2715]	89	!!
[503]	90	INTEGER, INTENT(in) :: kt ! ocean time-step index
[455]	91	!!----------------------------------------------------------------------
[2528]	92	!
[2715]	93	IF( wrk_in_use(3, 1,2) ) THEN
	94	CALL ctl_stop('dyn_hpg: requested workspace arrays are unavailable') ; RETURN
	95	ENDIF
	96	!
[503]	97	IF( l_trddyn ) THEN ! Temporary saving of ua and va trends (l_trddyn)
[455]	98	ztrdu(:,:,:) = ua(:,:,:)
	99	ztrdv(:,:,:) = va(:,:,:)
	100	ENDIF
[2528]	101	!
[455]	102	SELECT CASE ( nhpg ) ! Hydrastatic pressure gradient computation
[503]	103	CASE ( 0 ) ; CALL hpg_zco ( kt ) ! z-coordinate
	104	CASE ( 1 ) ; CALL hpg_zps ( kt ) ! z-coordinate plus partial steps (interpolation)
	105	CASE ( 2 ) ; CALL hpg_sco ( kt ) ! s-coordinate (standard jacobian formulation)
	106	CASE ( 3 ) ; CALL hpg_hel ( kt ) ! s-coordinate (helsinki modification)
	107	CASE ( 4 ) ; CALL hpg_wdj ( kt ) ! s-coordinate (weighted density jacobian)
	108	CASE ( 5 ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial)
	109	CASE ( 6 ) ; CALL hpg_rot ( kt ) ! s-coordinate (ROTated axes scheme)
[455]	110	END SELECT
[2528]	111	!
[503]	112	IF( l_trddyn ) THEN ! save the hydrostatic pressure gradient trends for momentum trend diagnostics
[455]	113	ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:)
	114	ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:)
[503]	115	CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )
[455]	116	ENDIF
[503]	117	!
	118	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' hpg - Ua: ', mask1=umask, &
	119	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
	120	!
[2715]	121	IF( wrk_not_released(3, 1,2) ) CALL ctl_stop('dyn_hpg: failed to release workspace arrays')
	122	!
[455]	123	END SUBROUTINE dyn_hpg
	124
	125
[2528]	126	SUBROUTINE dyn_hpg_init
[455]	127	!!----------------------------------------------------------------------
[2528]	128	!! * ROUTINE dyn_hpg_init *
[455]	129	!!
	130	!! ** Purpose : initializations for the hydrostatic pressure gradient
	131	!! computation and consistency control
	132	!!
[1601]	133	!! ** Action : Read the namelist namdyn_hpg and check the consistency
[455]	134	!! with the type of vertical coordinate used (zco, zps, sco)
	135	!!----------------------------------------------------------------------
	136	INTEGER :: ioptio = 0 ! temporary integer
[1601]	137	!!
[2528]	138	NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_hel, &
	139	& ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, rn_gamma , ln_dynhpg_imp
[455]	140	!!----------------------------------------------------------------------
[2528]	141	!
	142	REWIND( numnam ) ! Read Namelist namdyn_hpg
	143	READ ( numnam, namdyn_hpg )
	144	!
	145	IF(lwp) THEN ! Control print
[455]	146	WRITE(numout,*)
[2528]	147	WRITE(numout,*) 'dyn_hpg_init : hydrostatic pressure gradient initialisation'
	148	WRITE(numout,*) '~~~~~~~~~~~~'
[1601]	149	WRITE(numout,*) ' Namelist namdyn_hpg : choice of hpg scheme'
	150	WRITE(numout,*) ' z-coord. - full steps ln_hpg_zco = ', ln_hpg_zco
	151	WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps
	152	WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco
	153	WRITE(numout,*) ' s-coord. (helsinki modification) ln_hpg_hel = ', ln_hpg_hel
	154	WRITE(numout,*) ' s-coord. (weighted density jacobian) ln_hpg_wdj = ', ln_hpg_wdj
	155	WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc
	156	WRITE(numout,*) ' s-coord. (ROTated axes scheme) ln_hpg_rot = ', ln_hpg_rot
	157	WRITE(numout,*) ' weighting coeff. (wdj scheme) rn_gamma = ', rn_gamma
	158	WRITE(numout,*) ' time stepping: centered (F) or semi-implicit (T) ln_dynhpg_imp = ', ln_dynhpg_imp
[455]	159	ENDIF
[2528]	160	!
	161	IF( lk_vvl .AND. .NOT. ln_hpg_sco ) &
[2715]	162	& CALL ctl_stop('dyn_hpg_init : variable volume key_vvl require the standard jacobian formulation hpg_sco')
[2528]	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
	168	IF( ln_hpg_hel ) nhpg = 3
	169	IF( ln_hpg_wdj ) nhpg = 4
	170	IF( ln_hpg_djc ) nhpg = 5
	171	IF( ln_hpg_rot ) nhpg = 6
[2528]	172	!
[503]	173	! ! Consitency check
[455]	174	ioptio = 0
	175	IF( ln_hpg_zco ) ioptio = ioptio + 1
	176	IF( ln_hpg_zps ) ioptio = ioptio + 1
	177	IF( ln_hpg_sco ) ioptio = ioptio + 1
	178	IF( ln_hpg_hel ) ioptio = ioptio + 1
	179	IF( ln_hpg_wdj ) ioptio = ioptio + 1
	180	IF( ln_hpg_djc ) ioptio = ioptio + 1
	181	IF( ln_hpg_rot ) ioptio = ioptio + 1
[2715]	182	IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' )
[503]	183	!
[2528]	184	END SUBROUTINE dyn_hpg_init
[455]	185
	186
	187	SUBROUTINE hpg_zco( kt )
	188	!!---------------------------------------------------------------------
	189	!! * ROUTINE hpg_zco *
	190	!!
	191	!! ** Method : z-coordinate case, levels are horizontal surfaces.
	192	!! The now hydrostatic pressure gradient at a given level, jk,
	193	!! is computed by taking the vertical integral of the in-situ
	194	!! density gradient along the model level from the suface to that
	195	!! level: zhpi = grav .....
	196	!! zhpj = grav .....
[3]	197	!! add it to the general momentum trend (ua,va).
[455]	198	!! ua = ua - 1/e1u * zhpi
	199	!! va = va - 1/e2v * zhpj
	200	!!
[3]	201	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
[503]	202	!!----------------------------------------------------------------------
[2715]	203	USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
[3211]	204	!! DCSE_NEMO: need additional directives for renamed module variables
	205	!FTRANS zhpi :I :I :z
	206	!FTRANS zhpj :I :I :z
	207
[3]	208	!!
[503]	209	INTEGER, INTENT(in) :: kt ! ocean time-step index
	210	!!
	211	INTEGER :: ji, jj, jk ! dummy loop indices
	212	REAL(wp) :: zcoef0, zcoef1 ! temporary scalars
[3]	213	!!----------------------------------------------------------------------
[455]	214
[3]	215	IF( kt == nit000 ) THEN
	216	IF(lwp) WRITE(numout,*)
[455]	217	IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend'
	218	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case '
[3]	219	ENDIF
[455]	220
[2715]	221	zcoef0 = - grav * 0.5_wp ! Local constant initialization
[3]	222
[455]	223	! Surface value
[3]	224	DO jj = 2, jpjm1
	225	DO ji = fs_2, fs_jpim1 ! vector opt.
[455]	226	zcoef1 = zcoef0 * fse3w(ji,jj,1)
	227	! hydrostatic pressure gradient
	228	zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
	229	zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
[3]	230	! add to the general momentum trend
[455]	231	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
	232	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
	233	END DO
	234	END DO
[503]	235	!
[455]	236	! interior value (2=<jk=<jpkm1)
[3211]	237	#if defined key_z_first
	238	DO jj = 2, jpjm1
	239	DO ji = 2, jpim1
	240	DO jk = 2, jpkm1
	241	#else
[3]	242	DO jk = 2, jpkm1
[455]	243	DO jj = 2, jpjm1
[3]	244	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	245	#endif
[455]	246	zcoef1 = zcoef0 * fse3w(ji,jj,jk)
	247	! hydrostatic pressure gradient
	248	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &
	249	& + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) &
	250	& - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) / e1u(ji,jj)
	251
	252	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &
	253	& + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) &
	254	& - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) / e2v(ji,jj)
[3]	255	! add to the general momentum trend
[455]	256	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
	257	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
[3]	258	END DO
	259	END DO
	260	END DO
[503]	261	!
[455]	262	END SUBROUTINE hpg_zco
[216]	263
[3]	264
[455]	265	SUBROUTINE hpg_zps( kt )
[3]	266	!!---------------------------------------------------------------------
[455]	267	!! * ROUTINE hpg_zps *
[3]	268	!!
[455]	269	!! ** Method : z-coordinate plus partial steps case. blahblah...
	270	!!
[3]	271	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
[455]	272	!!----------------------------------------------------------------------
[2715]	273	USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
[3211]	274	!! DCSE_NEMO: need additional directives for renamed module variables
	275	!FTRANS zhpi :I :I :z
	276	!FTRANS zhpj :I :I :z
[503]	277	!!
	278	INTEGER, INTENT(in) :: kt ! ocean time-step index
	279	!!
	280	INTEGER :: ji, jj, jk ! dummy loop indices
	281	INTEGER :: iku, ikv ! temporary integers
	282	REAL(wp) :: zcoef0, zcoef1, zcoef2, zcoef3 ! temporary scalars
[3]	283	!!----------------------------------------------------------------------
	284
	285	IF( kt == nit000 ) THEN
	286	IF(lwp) WRITE(numout,*)
[455]	287	IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend'
[503]	288	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate with partial steps - vector optimization'
[3]	289	ENDIF
	290
[503]	291	! Local constant initialization
[2528]	292	zcoef0 = - grav * 0.5_wp
[3]	293
[2528]	294	! Surface value (also valid in partial step case)
[3]	295	DO jj = 2, jpjm1
	296	DO ji = fs_2, fs_jpim1 ! vector opt.
[170]	297	zcoef1 = zcoef0 * fse3w(ji,jj,1)
[3]	298	! hydrostatic pressure gradient
[455]	299	zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) / e1u(ji,jj)
	300	zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj)
[3]	301	! add to the general momentum trend
	302	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
	303	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
	304	END DO
	305	END DO
	306
[503]	307	! interior value (2=<jk=<jpkm1)
[3211]	308	#if defined key_z_first
	309	DO jj = 2, jpjm1
	310	DO ji = 2, jpim1
	311	DO jk = 2, jpkm1
	312	#else
[3]	313	DO jk = 2, jpkm1
	314	DO jj = 2, jpjm1
	315	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	316	#endif
[170]	317	zcoef1 = zcoef0 * fse3w(ji,jj,jk)
[3]	318	! hydrostatic pressure gradient
	319	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &
[455]	320	& + zcoef1 * ( ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) &
	321	& - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) / e1u(ji,jj)
[3]	322
	323	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &
[455]	324	& + zcoef1 * ( ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) &
	325	& - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) / e2v(ji,jj)
[3]	326	! add to the general momentum trend
	327	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
	328	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
[455]	329	END DO
[3]	330	END DO
	331	END DO
	332
[2528]	333	! partial steps correction at the last level (use gru & grv computed in zpshde.F90)
[3]	334	# if defined key_vectopt_loop
	335	jj = 1
	336	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
	337	# else
	338	DO jj = 2, jpjm1
	339	DO ji = 2, jpim1
	340	# endif
[2528]	341	iku = mbku(ji,jj)
	342	ikv = mbkv(ji,jj)
[3]	343	zcoef2 = zcoef0 * MIN( fse3w(ji,jj,iku), fse3w(ji+1,jj ,iku) )
	344	zcoef3 = zcoef0 * MIN( fse3w(ji,jj,ikv), fse3w(ji ,jj+1,ikv) )
[2528]	345	IF( iku > 1 ) THEN ! on i-direction (level 2 or more)
	346	ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value
	347	zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one
	348	& + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj)
	349	ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend
[3]	350	ENDIF
[2528]	351	IF( ikv > 1 ) THEN ! on j-direction (level 2 or more)
	352	va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value
	353	zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one
	354	& + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj)
	355	va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend
[3]	356	ENDIF
	357	# if ! defined key_vectopt_loop
	358	END DO
	359	# endif
	360	END DO
[503]	361	!
[455]	362	END SUBROUTINE hpg_zps
[216]	363
[3]	364
[455]	365	SUBROUTINE hpg_sco( kt )
[3]	366	!!---------------------------------------------------------------------
[455]	367	!! * ROUTINE hpg_sco *
[3]	368	!!
[455]	369	!! ** Method : s-coordinate case. Jacobian scheme.
	370	!! The now hydrostatic pressure gradient at a given level, jk,
	371	!! is computed by taking the vertical integral of the in-situ
[3]	372	!! density gradient along the model level from the suface to that
[455]	373	!! level. s-coordinates (ln_sco): a corrective term is added
	374	!! to the horizontal pressure gradient :
	375	!! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ]
	376	!! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ]
[3]	377	!! add it to the general momentum trend (ua,va).
[455]	378	!! ua = ua - 1/e1u * zhpi
	379	!! va = va - 1/e2v * zhpj
[3]	380	!!
	381	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
[503]	382	!!----------------------------------------------------------------------
[2715]	383	USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
[3211]	384	!! DCSE_NEMO: need additional directives for renamed module variables
	385	!FTRANS zhpi :I :I :z
	386	!FTRANS zhpj :I :I :z
[3]	387	!!
[503]	388	INTEGER, INTENT(in) :: kt ! ocean time-step index
	389	!!
[592]	390	INTEGER :: ji, jj, jk ! dummy loop indices
	391	REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars
[3]	392	!!----------------------------------------------------------------------
	393
	394	IF( kt == nit000 ) THEN
	395	IF(lwp) WRITE(numout,*)
[455]	396	IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend'
	397	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used'
[3]	398	ENDIF
	399
[503]	400	! Local constant initialization
[2528]	401	zcoef0 = - grav * 0.5_wp
[592]	402	! To use density and not density anomaly
[2528]	403	IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume
	404	ELSE ; znad = 0._wp ! Fixed volume
[592]	405	ENDIF
[455]	406
[503]	407	! Surface value
[455]	408	DO jj = 2, jpjm1
	409	DO ji = fs_2, fs_jpim1 ! vector opt.
	410	! hydrostatic pressure gradient along s-surfaces
[592]	411	zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) &
	412	& - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) )
	413	zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) &
	414	& - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) )
[455]	415	! s-coordinate pressure gradient correction
[2528]	416	zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) &
[455]	417	& * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
[2528]	418	zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) &
[455]	419	& * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
	420	! add to the general momentum trend
	421	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
	422	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
	423	END DO
	424	END DO
	425
[503]	426	! interior value (2=<jk=<jpkm1)
[3211]	427	#if defined key_z_first
	428	DO jj = 2, jpjm1
	429	DO ji = 2, jpim1
	430	DO jk = 2, jpkm1
	431	#else
[455]	432	DO jk = 2, jpkm1
	433	DO jj = 2, jpjm1
	434	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	435	#endif
[455]	436	! hydrostatic pressure gradient along s-surfaces
	437	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) &
[592]	438	& * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) &
	439	& - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) )
[455]	440	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) &
[592]	441	& * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) &
	442	& - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) )
[455]	443	! s-coordinate pressure gradient correction
[2528]	444	zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) &
[455]	445	& * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj)
[2528]	446	zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) &
[455]	447	& * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj)
	448	! add to the general momentum trend
	449	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
	450	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
	451	END DO
	452	END DO
	453	END DO
[503]	454	!
[455]	455	END SUBROUTINE hpg_sco
	456
	457
	458	SUBROUTINE hpg_hel( kt )
	459	!!---------------------------------------------------------------------
	460	!! * ROUTINE hpg_hel *
	461	!!
	462	!! ** Method : s-coordinate case.
	463	!! The now hydrostatic pressure gradient at a given level
	464	!! jk is computed by taking the vertical integral of the in-situ
	465	!! density gradient along the model level from the suface to that
	466	!! level. s-coordinates (ln_sco): a corrective term is added
	467	!! to the horizontal pressure gradient :
	468	!! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ]
	469	!! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ]
	470	!! add it to the general momentum trend (ua,va).
	471	!! ua = ua - 1/e1u * zhpi
	472	!! va = va - 1/e2v * zhpj
	473	!!
	474	!! ** Action : - Update (ua,va) with the now hydrastatic pressure trend
	475	!! - Save the trend (l_trddyn=T)
[503]	476	!!----------------------------------------------------------------------
[2715]	477	USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
[3211]	478	!! DCSE_NEMO: need additional directives for renamed module variables
	479	!FTRANS zhpi :I :I :z
	480	!FTRANS zhpj :I :I :z
[455]	481	!!
[503]	482	INTEGER, INTENT(in) :: kt ! ocean time-step index
	483	!!
	484	INTEGER :: ji, jj, jk ! dummy loop indices
	485	REAL(wp) :: zcoef0, zuap, zvap ! temporary scalars
[455]	486	!!----------------------------------------------------------------------
	487
	488	IF( kt == nit000 ) THEN
	489	IF(lwp) WRITE(numout,*)
	490	IF(lwp) WRITE(numout,*) 'dyn:hpg_hel : hydrostatic pressure gradient trend'
	491	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, helsinki modified scheme'
[216]	492	ENDIF
	493
[503]	494	! Local constant initialization
[2528]	495	zcoef0 = - grav * 0.5_wp
[455]	496
[503]	497	! Surface value
[3]	498	DO jj = 2, jpjm1
	499	DO ji = fs_2, fs_jpim1 ! vector opt.
[455]	500	! hydrostatic pressure gradient along s-surfaces
	501	zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj ,1) * rhd(ji+1,jj ,1) &
	502	& - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) )
	503	zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji ,jj+1,1) * rhd(ji ,jj+1,1) &
	504	& - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) )
	505	! s-coordinate pressure gradient correction
	506	zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) &
	507	& * ( fsdept(ji+1,jj,1) - fsdept(ji,jj,1) ) / e1u(ji,jj)
	508	zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) &
	509	& * ( fsdept(ji,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj)
[3]	510	! add to the general momentum trend
[455]	511	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
	512	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
[3]	513	END DO
	514	END DO
[503]	515	!
	516	! interior value (2=<jk=<jpkm1)
[3211]	517	#if defined key_z_first
	518	DO jj = 2, jpjm1
	519	DO ji = 2, jpim1
	520	DO jk = 2, jpkm1
	521	#else
[3]	522	DO jk = 2, jpkm1
	523	DO jj = 2, jpjm1
	524	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	525	#endif
[455]	526	! hydrostatic pressure gradient along s-surfaces
	527	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &
	528	& + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk ) * rhd(ji+1,jj,jk) &
	529	& -fse3t(ji ,jj,jk ) * rhd(ji ,jj,jk) ) &
	530	& + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1) &
	531	& -fse3t(ji ,jj,jk-1) * rhd(ji ,jj,jk-1) )
	532	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &
[2528]	533	& + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk) &
[455]	534	& -fse3t(ji,jj ,jk ) * rhd(ji,jj, jk) ) &
[2528]	535	& + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) &
[455]	536	& -fse3t(ji,jj ,jk-1) * rhd(ji,jj, jk-1) )
	537	! s-coordinate pressure gradient correction
	538	zuap = - zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) &
	539	& * ( fsdept(ji+1,jj,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)
	540	zvap = - zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) &
	541	& * ( fsdept(ji,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)
	542	! add to the general momentum trend
	543	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap
	544	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap
	545	END DO
	546	END DO
	547	END DO
[503]	548	!
[455]	549	END SUBROUTINE hpg_hel
	550
	551
	552	SUBROUTINE hpg_wdj( kt )
	553	!!---------------------------------------------------------------------
	554	!! * ROUTINE hpg_wdj *
	555	!!
	556	!! ** Method : Weighted Density Jacobian (wdj) scheme (song 1998)
[1601]	557	!! The weighting coefficients from the namelist parameter rn_gamma
	558	!! (alpha=0.5-rn_gamma ; beta=1-alpha=0.5+rn_gamma
[455]	559	!!
	560	!! Reference : Song, Mon. Wea. Rev., 126, 3213-3230, 1998.
[503]	561	!!----------------------------------------------------------------------
[2715]	562	USE oce, ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
[3211]	563	!! DCSE_NEMO: need additional directives for renamed module variables
	564	!FTRANS zhpi :I :I :z
	565	!FTRANS zhpj :I :I :z
[455]	566	!!
[503]	567	INTEGER, INTENT(in) :: kt ! ocean time-step index
	568	!!
	569	INTEGER :: ji, jj, jk ! dummy loop indices
	570	REAL(wp) :: zcoef0, zuap, zvap ! temporary scalars
	571	REAL(wp) :: zalph , zbeta ! " "
[455]	572	!!----------------------------------------------------------------------
	573
	574	IF( kt == nit000 ) THEN
	575	IF(lwp) WRITE(numout,*)
	576	IF(lwp) WRITE(numout,*) 'dyn:hpg_wdj : hydrostatic pressure gradient trend'
	577	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ Weighted Density Jacobian'
	578	ENDIF
	579
	580	! Local constant initialization
[2528]	581	zcoef0 = - grav * 0.5_wp
	582	zalph = 0.5_wp - rn_gamma ! weighting coefficients (alpha=0.5-rn_gamma
	583	zbeta = 0.5_wp + rn_gamma ! (beta =1-alpha=0.5+rn_gamma
[455]	584
	585	! Surface value (no ponderation)
	586	DO jj = 2, jpjm1
	587	DO ji = fs_2, fs_jpim1 ! vector opt.
	588	! hydrostatic pressure gradient along s-surfaces
	589	zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * rhd(ji+1,jj ,1) &
	590	& - fse3w(ji ,jj ,1) * rhd(ji ,jj ,1) )
	591	zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * rhd(ji ,jj+1,1) &
	592	& - fse3w(ji ,jj ,1) * rhd(ji, jj ,1) )
	593	! s-coordinate pressure gradient correction
	594	zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) &
	595	& * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj)
	596	zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) &
	597	& * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj)
	598	! add to the general momentum trend
	599	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap
	600	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap
	601	END DO
	602	END DO
	603
	604	! Interior value (2=<jk=<jpkm1) (weighted with zalph & zbeta)
[3211]	605	#if defined key_z_first
	606	DO jj = 2, jpjm1
	607	DO ji = 2, jpim1
	608	DO jk = 2, jpkm1
	609	#else
[455]	610	DO jk = 2, jpkm1
	611	DO jj = 2, jpjm1
	612	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	613	#endif
[455]	614	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) &
	615	& * ( ( fsde3w(ji+1,jj,jk ) + fsde3w(ji,jj,jk ) &
	616	& - fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) ) &
	617	& * ( zalph * ( rhd (ji+1,jj,jk-1) - rhd (ji,jj,jk-1) ) &
	618	& + zbeta * ( rhd (ji+1,jj,jk ) - rhd (ji,jj,jk ) ) ) &
	619	& - ( rhd (ji+1,jj,jk ) + rhd (ji,jj,jk ) &
	620	& - rhd (ji+1,jj,jk-1) - rhd (ji,jj,jk-1) ) &
	621	& * ( zalph * ( fsde3w(ji+1,jj,jk-1) - fsde3w(ji,jj,jk-1) ) &
	622	& + zbeta * ( fsde3w(ji+1,jj,jk ) - fsde3w(ji,jj,jk ) ) ) )
	623	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) &
	624	& * ( ( fsde3w(ji,jj+1,jk ) + fsde3w(ji,jj,jk ) &
	625	& - fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) ) &
	626	& * ( zalph * ( rhd (ji,jj+1,jk-1) - rhd (ji,jj,jk-1) ) &
	627	& + zbeta * ( rhd (ji,jj+1,jk ) - rhd (ji,jj,jk ) ) ) &
	628	& - ( rhd (ji,jj+1,jk ) + rhd (ji,jj,jk ) &
	629	& - rhd (ji,jj+1,jk-1) - rhd (ji,jj,jk-1) ) &
	630	& * ( zalph * ( fsde3w(ji,jj+1,jk-1) - fsde3w(ji,jj,jk-1) ) &
	631	& + zbeta * ( fsde3w(ji,jj+1,jk ) - fsde3w(ji,jj,jk ) ) ) )
[3]	632	! add to the general momentum trend
	633	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
	634	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
[455]	635	END DO
[3]	636	END DO
	637	END DO
[503]	638	!
[455]	639	END SUBROUTINE hpg_wdj
[216]	640
[455]	641
	642	SUBROUTINE hpg_djc( kt )
	643	!!---------------------------------------------------------------------
	644	!! * ROUTINE hpg_djc *
	645	!!
	646	!! ** Method : Density Jacobian with Cubic polynomial scheme
	647	!!
[503]	648	!! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003
[455]	649	!!----------------------------------------------------------------------
[2715]	650	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	651	USE oce , ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
	652	USE wrk_nemo, ONLY: drhox => wrk_3d_1 , dzx => wrk_3d_2
	653	USE wrk_nemo, ONLY: drhou => wrk_3d_3 , dzu => wrk_3d_4 , rho_i => wrk_3d_5
	654	USE wrk_nemo, ONLY: drhoy => wrk_3d_6 , dzy => wrk_3d_7
	655	USE wrk_nemo, ONLY: drhov => wrk_3d_8 , dzv => wrk_3d_9 , rho_j => wrk_3d_10
	656	USE wrk_nemo, ONLY: drhoz => wrk_3d_11 , dzz => wrk_3d_12
	657	USE wrk_nemo, ONLY: drhow => wrk_3d_13 , dzw => wrk_3d_14
	658	USE wrk_nemo, ONLY: rho_k => wrk_3d_15
[3211]	659	!! DCSE_NEMO: need additional directives for renamed module variables
	660	!FTRANS zhpi :I :I :z
	661	!FTRANS zhpj :I :I :z
	662	!FTRANS drhox :I :I :z
	663	!FTRANS dzx :I :I :z
	664	!FTRANS drhou :I :I :z
	665	!FTRANS dzu :I :I :z
	666	!FTRANS rho_i :I :I :z
	667	!FTRANS drhoy :I :I :z
	668	!FTRANS dzy :I :I :z
	669	!FTRANS drhov :I :I :z
	670	!FTRANS dzv :I :I :z
	671	!FTRANS rho_j :I :I :z
	672	!FTRANS drhoz :I :I :z
	673	!FTRANS dzz :I :I :z
	674	!FTRANS drhow :I :I :z
	675	!FTRANS dzw :I :I :z
	676	!FTRANS rho_k :I :I :z
[503]	677	!!
	678	INTEGER, INTENT(in) :: kt ! ocean time-step index
	679	!!
	680	INTEGER :: ji, jj, jk ! dummy loop indices
	681	REAL(wp) :: zcoef0, zep, cffw ! temporary scalars
	682	REAL(wp) :: z1_10, cffu, cffx ! " "
	683	REAL(wp) :: z1_12, cffv, cffy ! " "
[455]	684	!!----------------------------------------------------------------------
	685
[2715]	686	IF( wrk_in_use(3, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15) ) THEN
	687	CALL ctl_stop('dyn:hpg_djc: requested workspace arrays unavailable') ; RETURN
	688	ENDIF
	689
[455]	690	IF( kt == nit000 ) THEN
	691	IF(lwp) WRITE(numout,*)
	692	IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend'
	693	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme'
[216]	694	ENDIF
	695
[503]	696	! Local constant initialization
[2528]	697	zcoef0 = - grav * 0.5_wp
	698	z1_10 = 1._wp / 10._wp
	699	z1_12 = 1._wp / 12._wp
[455]	700
	701	!----------------------------------------------------------------------------------------
	702	! compute and store in provisional arrays elementary vertical and horizontal differences
	703	!----------------------------------------------------------------------------------------
	704
	705	!!bug gm Not a true bug, but... dzz=e3w for dzx, dzy verify what it is really
	706
[3211]	707	#if defined key_z_first
	708	DO jj = 2, jpjm1
	709	DO ji = 2, jpim1
	710	DO jk = 2, jpkm1
	711	#else
[455]	712	DO jk = 2, jpkm1
	713	DO jj = 2, jpjm1
	714	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	715	#endif
[455]	716	drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1)
	717	dzz (ji,jj,jk) = fsde3w(ji ,jj ,jk) - fsde3w(ji,jj,jk-1)
	718	drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk )
	719	dzx (ji,jj,jk) = fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk )
	720	drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk )
	721	dzy (ji,jj,jk) = fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk )
	722	END DO
	723	END DO
	724	END DO
	725
	726	!-------------------------------------------------------------------------
	727	! compute harmonic averages using eq. 5.18
	728	!-------------------------------------------------------------------------
	729	zep = 1.e-15
	730
[503]	731	!!bug gm drhoz not defined at level 1 and used (jk-1 with jk=2)
	732	!!bug gm idem for drhox, drhoy et ji=jpi and jj=jpj
[455]	733
[3211]	734	#if defined key_z_first
	735	DO jj = 2, jpjm1
	736	DO ji = 2, jpim1
	737	DO jk = 2, jpkm1
	738	#else
[455]	739	DO jk = 2, jpkm1
	740	DO jj = 2, jpjm1
	741	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	742	#endif
[2528]	743	cffw = 2._wp * drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1)
[455]	744
[2528]	745	cffu = 2._wp * drhox(ji+1,jj ,jk) * drhox(ji,jj,jk )
	746	cffx = 2._wp * dzx (ji+1,jj ,jk) * dzx (ji,jj,jk )
[455]	747
[2528]	748	cffv = 2._wp * drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk )
	749	cffy = 2._wp * dzy (ji ,jj+1,jk) * dzy (ji,jj,jk )
[455]	750
	751	IF( cffw > zep) THEN
[2528]	752	drhow(ji,jj,jk) = 2._wp * drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1) &
	753	& / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) )
[455]	754	ELSE
[2528]	755	drhow(ji,jj,jk) = 0._wp
[455]	756	ENDIF
	757
[2528]	758	dzw(ji,jj,jk) = 2._wp * dzz(ji,jj,jk) * dzz(ji,jj,jk-1) &
	759	& / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) )
[455]	760
	761	IF( cffu > zep ) THEN
[2528]	762	drhou(ji,jj,jk) = 2._wp * drhox(ji+1,jj,jk) * drhox(ji,jj,jk) &
	763	& / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) )
[455]	764	ELSE
[2528]	765	drhou(ji,jj,jk ) = 0._wp
[455]	766	ENDIF
	767
	768	IF( cffx > zep ) THEN
[2528]	769	dzu(ji,jj,jk) = 2._wp * dzx(ji+1,jj,jk) * dzx(ji,jj,jk) &
	770	& / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) )
[455]	771	ELSE
[2528]	772	dzu(ji,jj,jk) = 0._wp
[455]	773	ENDIF
	774
	775	IF( cffv > zep ) THEN
[2528]	776	drhov(ji,jj,jk) = 2._wp * drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk) &
	777	& / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) )
[455]	778	ELSE
[2528]	779	drhov(ji,jj,jk) = 0._wp
[455]	780	ENDIF
	781
	782	IF( cffy > zep ) THEN
[2528]	783	dzv(ji,jj,jk) = 2._wp * dzy(ji,jj+1,jk) * dzy(ji,jj,jk) &
	784	& / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) )
[455]	785	ELSE
[2528]	786	dzv(ji,jj,jk) = 0._wp
[455]	787	ENDIF
	788
	789	END DO
	790	END DO
	791	END DO
	792
	793	!----------------------------------------------------------------------------------
	794	! apply boundary conditions at top and bottom using 5.36-5.37
	795	!----------------------------------------------------------------------------------
[2528]	796	drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:, 1 ) ) - 0.5_wp * drhow(:,:, 2 )
	797	drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:, 1 ) ) - 0.5_wp * drhou(:,:, 2 )
	798	drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:, 1 ) ) - 0.5_wp * drhov(:,:, 2 )
[455]	799
[2528]	800	drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1)
	801	drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1)
	802	drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1)
[455]	803
	804
	805	!--------------------------------------------------------------
	806	! Upper half of top-most grid box, compute and store
	807	!-------------------------------------------------------------
	808
	809	!!bug gm : e3w-de3w = 0.5*e3w .... and de3w(2)-de3w(1)=e3w(2) .... to be verified
	810	! true if de3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be
	811
	812	DO jj = 2, jpjm1
	813	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	814	rho_k(ji,jj,1) = -grav * ( fse3w(ji,jj,1) - fsde3w(ji,jj,1) ) &
	815	& * ( rhd(ji,jj,1) &
	816	& + 0.5_wp * ( rhd(ji,jj,2) - rhd(ji,jj,1) ) &
	817	& * ( fse3w (ji,jj,1) - fsde3w(ji,jj,1) ) &
	818	& / ( fsde3w(ji,jj,2) - fsde3w(ji,jj,1) ) )
[455]	819	END DO
	820	END DO
	821
	822	!!bug gm : here also, simplification is possible
	823	!!bug gm : optimisation: 1/10 and 1/12 the division should be done before the loop
	824
[3211]	825	#if defined key_z_first
	826	DO jj = 2, jpjm1
	827	DO ji = 2, jpim1
	828	DO jk = 2, jpkm1
	829	#else
[455]	830	DO jk = 2, jpkm1
	831	DO jj = 2, jpjm1
	832	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	833	#endif
[455]	834
	835	rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) &
	836	& * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) ) &
	837	& - grav * z1_10 * ( &
	838	& ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) &
	839	& * ( fsde3w(ji,jj,jk) - fsde3w(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) &
	840	& - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) &
	841	& * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) &
	842	& )
	843
	844	rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) &
	845	& * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) ) &
	846	& - grav* z1_10 * ( &
	847	& ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) &
	848	& * ( fsde3w(ji+1,jj,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) &
	849	& - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) &
	850	& * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) &
	851	& )
	852
	853	rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) &
	854	& * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) ) &
	855	& - grav* z1_10 * ( &
	856	& ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) &
	857	& * ( fsde3w(ji,jj+1,jk) - fsde3w(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) &
	858	& - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) &
	859	& * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) &
	860	& )
	861
	862	END DO
	863	END DO
	864	END DO
	865	CALL lbc_lnk(rho_k,'W',1.)
	866	CALL lbc_lnk(rho_i,'U',1.)
	867	CALL lbc_lnk(rho_j,'V',1.)
	868
	869
	870	! ---------------
	871	! Surface value
	872	! ---------------
	873	DO jj = 2, jpjm1
	874	DO ji = fs_2, fs_jpim1 ! vector opt.
	875	zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj)
	876	zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj)
	877	! add to the general momentum trend
	878	ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1)
	879	va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1)
	880	END DO
	881	END DO
	882
	883	! ----------------
	884	! interior value (2=<jk=<jpkm1)
	885	! ----------------
[3211]	886	#if defined key_z_first
	887	DO jj = 2, jpjm1
	888	DO ji = 2, jpim1
	889	DO jk = 2, jpkm1
	890	#else
[455]	891	DO jk = 2, jpkm1
	892	DO jj = 2, jpjm1
	893	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	894	#endif
[455]	895	! hydrostatic pressure gradient along s-surfaces
	896	zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &
	897	& + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) ) &
	898	& - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) / e1u(ji,jj)
	899	zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &
	900	& + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) &
	901	& -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) / e2v(ji,jj)
	902	! add to the general momentum trend
	903	ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk)
	904	va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk)
	905	END DO
	906	END DO
	907	END DO
[503]	908	!
[2715]	909	IF( wrk_not_released(3, 1,2,3,4,5,6,7,8,9,10,11,12,13,14,15) ) &
	910	CALL ctl_stop('dyn:hpg_djc: failed to release workspace arrays')
	911	!
[455]	912	END SUBROUTINE hpg_djc
	913
	914
	915	SUBROUTINE hpg_rot( kt )
	916	!!---------------------------------------------------------------------
	917	!! * ROUTINE hpg_rot *
	918	!!
	919	!! ** Method : rotated axes scheme (Thiem and Berntsen 2005)
	920	!!
	921	!! Reference: Thiem & Berntsen, Ocean Modelling, In press, 2005.
	922	!!----------------------------------------------------------------------
[2715]	923	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	924	USE oce , ONLY: zhpi => ta , zhpj => sa ! (ta,sa) used as 3D workspace
	925	USE wrk_nemo, ONLY: zdistr => wrk_2d_1 , zsina => wrk_2d_2 , zcosa => wrk_2d_3
	926	USE wrk_nemo, ONLY: zhpiorg => wrk_3d_1 , zhpirot => wrk_3d_2
	927	USE wrk_nemo, ONLY: zhpitra => wrk_3d_3 , zhpine => wrk_3d_4
	928	USE wrk_nemo, ONLY: zhpjorg => wrk_3d_5 , zhpjrot => wrk_3d_6
	929	USE wrk_nemo, ONLY: zhpjtra => wrk_3d_7 , zhpjne => wrk_3d_8
[3211]	930	!! DCSE_NEMO: need additional directives for renamed module variables
	931	!FTRANS zhpi :I :I :z
	932	!FTRANS zhpj :I :I :z
	933	!FTRANS zhpiorg :I :I :z
	934	!FTRANS zhpirot :I :I :z
	935	!FTRANS zhpitra :I :I :z
	936	!FTRANS zhpine :I :I :z
	937	!FTRANS zhpjorg :I :I :z
	938	!FTRANS zhpjrot :I :I :z
	939	!FTRANS zhpjtra :I :I :z
	940	!FTRANS zhpjne :I :I :z
[503]	941	!!
	942	INTEGER, INTENT(in) :: kt ! ocean time-step index
	943	!!
	944	INTEGER :: ji, jj, jk ! dummy loop indices
	945	REAL(wp) :: zforg, zcoef0, zuap, zmskd1, zmskd1m ! temporary scalar
	946	REAL(wp) :: zfrot , zvap, zmskd2, zmskd2m ! " "
[455]	947	!!----------------------------------------------------------------------
	948
[2715]	949	IF( wrk_in_use(2, 1,2,3) .OR. &
	950	wrk_in_use(3, 1,2,3,4,5,6,7,8) ) THEN
	951	CALL ctl_stop('dyn:hpg_rot: requested workspace arrays unavailable') ; RETURN
	952	ENDIF
	953
[455]	954	IF( kt == nit000 ) THEN
	955	IF(lwp) WRITE(numout,*)
	956	IF(lwp) WRITE(numout,*) 'dyn:hpg_rot : hydrostatic pressure gradient trend'
	957	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, rotated axes scheme used'
[3]	958	ENDIF
	959
[455]	960	! -------------------------------
	961	! Local constant initialization
	962	! -------------------------------
[2528]	963	zcoef0 = - grav * 0.5_wp
	964	zforg = 0.95_wp
	965	zfrot = 1._wp - zforg
[3]	966
[455]	967	! inverse of the distance between 2 diagonal T-points (defined at F-point) (here zcoef0/distance)
	968	zdistr(:,:) = zcoef0 / SQRT( e1f(:,:)e1f(:,:) + e2f(:,:)e1f(:,:) )
[3]	969
[455]	970	! sinus and cosinus of diagonal angle at F-point
	971	zsina(:,:) = ATAN2( e2f(:,:), e1f(:,:) )
	972	zcosa(:,:) = COS( zsina(:,:) )
	973	zsina(:,:) = SIN( zsina(:,:) )
	974
	975	! ---------------
	976	! Surface value
	977	! ---------------
	978	! compute and add to the general trend the pressure gradients along the axes
	979	DO jj = 2, jpjm1
	980	DO ji = fs_2, fs_jpim1 ! vector opt.
	981	! hydrostatic pressure gradient along s-surfaces
	982	zhpiorg(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,1) * rhd(ji+1,jj,1) &
	983	& - fse3t(ji ,jj,1) * rhd(ji ,jj,1) )
	984	zhpjorg(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,1) * rhd(ji,jj+1,1) &
	985	& - fse3t(ji,jj ,1) * rhd(ji,jj ,1) )
	986	! s-coordinate pressure gradient correction
	987	zuap = -zcoef0 * ( rhd (ji+1,jj ,1) + rhd (ji,jj,1) ) &
	988	& * ( fsdept(ji+1,jj ,1) - fsdept(ji,jj,1) ) / e1u(ji,jj)
	989	zvap = -zcoef0 * ( rhd (ji ,jj+1,1) + rhd (ji,jj,1) ) &
	990	& * ( fsdept(ji ,jj+1,1) - fsdept(ji,jj,1) ) / e2v(ji,jj)
	991	! add to the general momentum trend
	992	ua(ji,jj,1) = ua(ji,jj,1) + zforg * ( zhpiorg(ji,jj,1) + zuap )
	993	va(ji,jj,1) = va(ji,jj,1) + zforg * ( zhpjorg(ji,jj,1) + zvap )
	994	END DO
	995	END DO
	996
	997	! compute the pressure gradients in the diagonal directions
	998	DO jj = 1, jpjm1
	999	DO ji = 1, fs_jpim1 ! vector opt.
[3211]	1000	#if defined key_z_first
	1001	zmskd1 = tmask_1(ji+1,jj+1) * tmask_1(ji ,jj) ! mask in the 1st diagnonal
	1002	zmskd2 = tmask_1(ji ,jj+1) * tmask_1(ji+1,jj) ! mask in the 2nd diagnonal
	1003	#else
[455]	1004	zmskd1 = tmask(ji+1,jj+1,1) * tmask(ji ,jj,1) ! mask in the 1st diagnonal
	1005	zmskd2 = tmask(ji ,jj+1,1) * tmask(ji+1,jj,1) ! mask in the 2nd diagnonal
[3211]	1006	#endif
[455]	1007	! hydrostatic pressure gradient along s-surfaces
	1008	zhpitra(ji,jj,1) = zdistr(ji,jj) * zmskd1 * ( fse3t(ji+1,jj+1,1) * rhd(ji+1,jj+1,1) &
	1009	& - fse3t(ji ,jj ,1) * rhd(ji ,jj ,1) )
	1010	zhpjtra(ji,jj,1) = zdistr(ji,jj) * zmskd2 * ( fse3t(ji ,jj+1,1) * rhd(ji ,jj+1,1) &
	1011	& - fse3t(ji+1,jj ,1) * rhd(ji+1,jj ,1) )
	1012	! s-coordinate pressure gradient correction
	1013	zuap = -zdistr(ji,jj) * zmskd1 * ( rhd (ji+1,jj+1,1) + rhd (ji ,jj,1) ) &
	1014	& * ( fsdept(ji+1,jj+1,1) - fsdept(ji ,jj,1) )
	1015	zvap = -zdistr(ji,jj) * zmskd2 * ( rhd (ji ,jj+1,1) + rhd (ji+1,jj,1) ) &
	1016	& * ( fsdept(ji ,jj+1,1) - fsdept(ji+1,jj,1) )
	1017	! back rotation
	1018	zhpine(ji,jj,1) = zcosa(ji,jj) * ( zhpitra(ji,jj,1) + zuap ) &
	1019	& - zsina(ji,jj) * ( zhpjtra(ji,jj,1) + zvap )
	1020	zhpjne(ji,jj,1) = zsina(ji,jj) * ( zhpitra(ji,jj,1) + zuap ) &
	1021	& + zcosa(ji,jj) * ( zhpjtra(ji,jj,1) + zvap )
	1022	END DO
	1023	END DO
	1024
	1025	! interpolate and add to the general trend the diagonal gradient
	1026	DO jj = 2, jpjm1
	1027	DO ji = fs_2, fs_jpim1 ! vector opt.
	1028	! averaging
	1029	zhpirot(ji,jj,1) = 0.5 * ( zhpine(ji,jj,1) + zhpine(ji ,jj-1,1) )
	1030	zhpjrot(ji,jj,1) = 0.5 * ( zhpjne(ji,jj,1) + zhpjne(ji-1,jj ,1) )
	1031	! add to the general momentum trend
	1032	ua(ji,jj,1) = ua(ji,jj,1) + zfrot * zhpirot(ji,jj,1)
	1033	va(ji,jj,1) = va(ji,jj,1) + zfrot * zhpjrot(ji,jj,1)
	1034	END DO
	1035	END DO
	1036
	1037	! -----------------
	1038	! 2. interior value (2=<jk=<jpkm1)
	1039	! -----------------
	1040	! compute and add to the general trend the pressure gradients along the axes
[3211]	1041	#if defined key_z_first
	1042	DO jj = 2, jpjm1
	1043	DO ji = 2, jpim1
	1044	DO jk = 2, jpkm1
	1045	#else
[455]	1046	DO jk = 2, jpkm1
	1047	DO jj = 2, jpjm1
	1048	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	1049	#endif
[455]	1050	! hydrostatic pressure gradient along s-surfaces
	1051	zhpiorg(ji,jj,jk) = zhpiorg(ji,jj,jk-1) &
	1052	& + zcoef0 / e1u(ji,jj) * ( fse3t(ji+1,jj,jk ) * rhd(ji+1,jj,jk ) &
	1053	& - fse3t(ji ,jj,jk ) * rhd(ji ,jj,jk ) &
	1054	& + fse3t(ji+1,jj,jk-1) * rhd(ji+1,jj,jk-1) &
	1055	& - fse3t(ji ,jj,jk-1) * rhd(ji ,jj,jk-1) )
	1056	zhpjorg(ji,jj,jk) = zhpjorg(ji,jj,jk-1) &
	1057	& + zcoef0 / e2v(ji,jj) * ( fse3t(ji,jj+1,jk ) * rhd(ji,jj+1,jk ) &
	1058	& - fse3t(ji,jj ,jk ) * rhd(ji,jj, jk ) &
	1059	& + fse3t(ji,jj+1,jk-1) * rhd(ji,jj+1,jk-1) &
	1060	& - fse3t(ji,jj ,jk-1) * rhd(ji,jj, jk-1) )
	1061	! s-coordinate pressure gradient correction
	1062	zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) ) &
	1063	& * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)
	1064	zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) ) &
	1065	& * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj)
	1066	! add to the general momentum trend
	1067	ua(ji,jj,jk) = ua(ji,jj,jk) + zforg*( zhpiorg(ji,jj,jk) + zuap )
	1068	va(ji,jj,jk) = va(ji,jj,jk) + zforg*( zhpjorg(ji,jj,jk) + zvap )
	1069	END DO
	1070	END DO
	1071	END DO
	1072
	1073	! compute the pressure gradients in the diagonal directions
[3211]	1074	#if defined key_z_first
	1075	DO jj = 1, jpjm1
	1076	DO ji = 1, jpim1
	1077	DO jk = 2, jpkm1
	1078	#else
[455]	1079	DO jk = 2, jpkm1
	1080	DO jj = 1, jpjm1
	1081	DO ji = 1, fs_jpim1 ! vector opt.
[3211]	1082	#endif
[455]	1083	zmskd1 = tmask(ji+1,jj+1,jk ) * tmask(ji ,jj,jk ) ! level jk mask in the 1st diagnonal
	1084	zmskd1m = tmask(ji+1,jj+1,jk-1) * tmask(ji ,jj,jk-1) ! level jk-1 " "
	1085	zmskd2 = tmask(ji ,jj+1,jk ) * tmask(ji+1,jj,jk ) ! level jk mask in the 2nd diagnonal
	1086	zmskd2m = tmask(ji ,jj+1,jk-1) * tmask(ji+1,jj,jk-1) ! level jk-1 " "
	1087	! hydrostatic pressure gradient along s-surfaces
	1088	zhpitra(ji,jj,jk) = zhpitra(ji,jj,jk-1) &
	1089	& + zdistr(ji,jj) * zmskd1 * ( fse3t(ji+1,jj+1,jk ) * rhd(ji+1,jj+1,jk) &
	1090	& -fse3t(ji ,jj ,jk ) * rhd(ji ,jj ,jk) ) &
	1091	& + zdistr(ji,jj) * zmskd1m * ( fse3t(ji+1,jj+1,jk-1) * rhd(ji+1,jj+1,jk-1) &
	1092	& -fse3t(ji ,jj ,jk-1) * rhd(ji ,jj ,jk-1) )
	1093	zhpjtra(ji,jj,jk) = zhpjtra(ji,jj,jk-1) &
	1094	& + zdistr(ji,jj) * zmskd2 * ( fse3t(ji ,jj+1,jk ) * rhd(ji ,jj+1,jk) &
	1095	& -fse3t(ji+1,jj ,jk ) * rhd(ji+1,jj, jk) ) &
	1096	& + zdistr(ji,jj) * zmskd2m * ( fse3t(ji ,jj+1,jk-1) * rhd(ji ,jj+1,jk-1) &
	1097	& -fse3t(ji+1,jj ,jk-1) * rhd(ji+1,jj, jk-1) )
	1098	! s-coordinate pressure gradient correction
	1099	zuap = - zdistr(ji,jj) * zmskd1 * ( rhd (ji+1,jj+1,jk) + rhd (ji ,jj,jk) ) &
	1100	& * ( fsdept(ji+1,jj+1,jk) - fsdept(ji ,jj,jk) )
	1101	zvap = - zdistr(ji,jj) * zmskd2 * ( rhd (ji ,jj+1,jk) + rhd (ji+1,jj,jk) ) &
	1102	& * ( fsdept(ji ,jj+1,jk) - fsdept(ji+1,jj,jk) )
	1103	! back rotation
	1104	zhpine(ji,jj,jk) = zcosa(ji,jj) * ( zhpitra(ji,jj,jk) + zuap ) &
	1105	& - zsina(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap )
	1106	zhpjne(ji,jj,jk) = zsina(ji,jj) * ( zhpitra(ji,jj,jk) + zuap ) &
	1107	& + zcosa(ji,jj) * ( zhpjtra(ji,jj,jk) + zvap )
	1108	END DO
	1109	END DO
	1110	END DO
	1111
	1112	! interpolate and add to the general trend
[3211]	1113	#if defined key_z_first
	1114	DO jj = 2, jpjm1
	1115	DO ji = 2, jpim1
	1116	DO jk = 2, jpkm1
	1117	#else
[455]	1118	DO jk = 2, jpkm1
	1119	DO jj = 2, jpjm1
	1120	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	1121	#endif
[455]	1122	! averaging
	1123	zhpirot(ji,jj,jk) = 0.5 * ( zhpine(ji,jj,jk) + zhpine(ji ,jj-1,jk) )
	1124	zhpjrot(ji,jj,jk) = 0.5 * ( zhpjne(ji,jj,jk) + zhpjne(ji-1,jj ,jk) )
	1125	! add to the general momentum trend
	1126	ua(ji,jj,jk) = ua(ji,jj,jk) + zfrot * zhpirot(ji,jj,jk)
	1127	va(ji,jj,jk) = va(ji,jj,jk) + zfrot * zhpjrot(ji,jj,jk)
	1128	END DO
	1129	END DO
	1130	END DO
[503]	1131	!
[2715]	1132	IF( wrk_not_released(2, 1,2,3) .OR. &
	1133	wrk_not_released(3, 1,2,3,4,5,6,7,8) ) CALL ctl_stop('dyn:hpg_rot: failed to release workspace arrays')
	1134	!
[455]	1135	END SUBROUTINE hpg_rot
	1136
[3]	1137	!!======================================================================
	1138	END MODULE dynhpg

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

Download in other formats: