Context Navigation

dynhpg.F90 @ 11101

Last change on this file since 11101 was 11101, checked in by frrh, 5 years ago

Merge changes from Met Office GMED ticket 450 to reduce unnecessary
text output from NEMO.
This output, which is typically not switchable, is rarely of interest
in normal (non-debugging) runs and simply redunantley consumes extra
file space.
Further, the presence of this text output has been shown to
significantly degrade performance of models which are run during
Met Office HPC RAID (disk) checks.
The new code introduces switches which are configurable via the
changes made in the associated Met Office MOCI ticket 399.

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 @ 11101

Download in other formats: