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

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 2621

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

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

Download in other formats: