Context Navigation

dynhpg.F90 @ 12396

Last change on this file since 12396 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Property svn:keywords set to Id

File size: 69.0 KB

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

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

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

Context Navigation

source: NEMO/trunk/src/OCE/DYN/dynhpg.F90 @ 12396

Download in other formats: