Context Navigation

dynhpg.F90 @ 3275

Last change on this file since 3275 was 2453, checked in by gm, 14 years ago
v3.3beta: #769 remove useless nn_dynhpg_rst namelist parameter
Property svn:keywords set to `Id`
File size: 51.9 KB

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 3275

Download in other formats: