Context Navigation

dynhpg.F90 @ 11872

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

Branch 2019/fix_sn_cfctl_ticket2328. See #2328. Replacement of ln_ctl and activation of full functionality with
sn_cfctl structure. These changes rename structure components l_mppout and l_mpptop as l_prtctl and l_prttrc
and introduce l_glochk to activate former ln_ctl code in stpctl.F90 to perform global location of min and max
checks. Also added is l_allon which can be used to activate all output (much like the former ln_ctl). If l_allon
is .false. then l_config decides whether or not the suboptions are used.

   sn_cfctl%l_glochk = .FALSE.    ! Range sanity checks are local (F) or global (T). Set T for debugging only
   sn_cfctl%l_allon  = .FALSE.    ! IF T activate all options. If F deactivate all unless l_config is T
   sn_cfctl%l_config = .TRUE.     ! IF .true. then control which reports are written with the remaining options

Note, these changes pass SETTE tests but all references to ln_ctl need to be removed from the sette scripts.

Property svn:keywords set to Id

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

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

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

Context Navigation

source: NEMO/branches/2019/fix_sn_cfctl_ticket2328/src/OCE/DYN/dynhpg.F90 @ 11872

Download in other formats: