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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

zpshde.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/zpshde.F90 @ 11057

Last change on this file since 11057 was 10954, checked in by acc, 5 years ago
2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert TRA modules and all knock on effects of these conversions. SETTE tested
Property svn:keywords set to `Id`
File size: 24.9 KB

Rev	Line
[3]	1	MODULE zpshde
[2528]	2	!!======================================================================
[3]	3	!! * MODULE zpshde *
[2528]	4	!! z-coordinate + partial step : Horizontal Derivative at ocean bottom level
	5	!!======================================================================
	6	!! History : OPA ! 2002-04 (A. Bozec) Original code
	7	!! NEMO 1.0 ! 2002-08 (G. Madec E. Durand) Optimization and Free form
	8	!! - ! 2004-03 (C. Ethe) adapted for passive tracers
	9	!! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA
[5120]	10	!! 3.6 ! 2014-11 (P. Mathiot) Add zps_hde_isf (needed to open a cavity)
[2528]	11	!!======================================================================
[457]	12
[3]	13	!!----------------------------------------------------------------------
	14	!! zps_hde : Horizontal DErivative of T, S and rd at the last
	15	!! ocean level (Z-coord. with Partial Steps)
	16	!!----------------------------------------------------------------------
[2528]	17	USE oce ! ocean: dynamics and tracers variables
	18	USE dom_oce ! domain: ocean variables
[3]	19	USE phycst ! physical constants
[2528]	20	USE eosbn2 ! ocean equation of state
[3]	21	USE in_out_manager ! I/O manager
	22	USE lbclnk ! lateral boundary conditions (or mpp link)
[2715]	23	USE lib_mpp ! MPP library
[3294]	24	USE timing ! Timing
[3]	25
	26	IMPLICIT NONE
	27	PRIVATE
	28
[5120]	29	PUBLIC zps_hde ! routine called by step.F90
	30	PUBLIC zps_hde_isf ! routine called by step.F90
[3]	31
	32	!! * Substitutions
	33	# include "vectopt_loop_substitute.h90"
	34	!!----------------------------------------------------------------------
[9598]	35	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[2528]	36	!! $Id$
[10068]	37	!! Software governed by the CeCILL license (see ./LICENSE)
[247]	38	!!----------------------------------------------------------------------
[3]	39	CONTAINS
	40
[10954]	41	SUBROUTINE zps_hde( kt, Kmm, kjpt, pta, pgtu, pgtv, &
[5120]	42	& prd, pgru, pgrv )
	43	!!----------------------------------------------------------------------
	44	!! * ROUTINE zps_hde *
	45	!!
	46	!! ** Purpose : Compute the horizontal derivative of T, S and rho
	47	!! at u- and v-points with a linear interpolation for z-coordinate
	48	!! with partial steps.
	49	!!
	50	!! ** Method : In z-coord with partial steps, scale factors on last
	51	!! levels are different for each grid point, so that T, S and rd
	52	!! points are not at the same depth as in z-coord. To have horizontal
	53	!! gradients again, we interpolate T and S at the good depth :
	54	!! Linear interpolation of T, S
	55	!! Computation of di(tb) and dj(tb) by vertical interpolation:
	56	!! di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~
	57	!! dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~
	58	!! This formulation computes the two cases:
	59	!! CASE 1 CASE 2
	60	!! k-1 ___ ___________ k-1 ___ ___________
	61	!! Ti T~ T~ Ti+1
	62	!! _____ _____
	63	!! k \| \|Ti+1 k Ti \| \|
	64	!! \| \|____ ____\| \|
	65	!! ___ \| \| \| ___ \| \| \|
	66	!!
	67	!! case 1-> e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
	68	!! t~ = t(i+1,j ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
	69	!! ( t~ = t(i ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1) )
	70	!! or
	71	!! case 2-> e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
	72	!! t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
	73	!! ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
	74	!! Idem for di(s) and dj(s)
	75	!!
	76	!! For rho, we call eos which will compute rd~(t~,s~) at the right
	77	!! depth zh from interpolated T and S for the different formulations
	78	!! of the equation of state (eos).
	79	!! Gradient formulation for rho :
	80	!! di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~
	81	!!
	82	!! ** Action : compute for top interfaces
	83	!! - pgtu, pgtv: horizontal gradient of tracer at u- & v-points
	84	!! - pgru, pgrv: horizontal gradient of rho (if present) at u- & v-points
	85	!!----------------------------------------------------------------------
	86	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[10954]	87	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[5120]	88	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	89	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields
	90	REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts
	91	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields
	92	REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom)
	93	!
[5836]	94	INTEGER :: ji, jj, jn ! Dummy loop indices
	95	INTEGER :: iku, ikv, ikum1, ikvm1 ! partial step level (ocean bottom level) at u- and v-points
	96	REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv ! local scalars
	97	REAL(wp), DIMENSION(jpi,jpj) :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos
	98	REAL(wp), DIMENSION(jpi,jpj,kjpt) :: zti, ztj !
[5120]	99	!!----------------------------------------------------------------------
	100	!
[9019]	101	IF( ln_timing ) CALL timing_start( 'zps_hde')
[5120]	102	!
[9019]	103	pgtu(:,:,:) = 0._wp ; zti (:,:,:) = 0._wp ; zhi (:,:) = 0._wp
	104	pgtv(:,:,:) = 0._wp ; ztj (:,:,:) = 0._wp ; zhj (:,:) = 0._wp
[5120]	105	!
	106	DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==!
	107	!
	108	DO jj = 1, jpjm1
	109	DO ji = 1, jpim1
	110	iku = mbku(ji,jj) ; ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points
	111	ikv = mbkv(ji,jj) ; ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1
[10954]	112	!!gm BUG ? when applied to before fields, e3w(:,:,:,Kbb) should be used....
	113	ze3wu = e3w(ji+1,jj ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
	114	ze3wv = e3w(ji ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
[5120]	115	!
	116	! i- direction
	117	IF( ze3wu >= 0._wp ) THEN ! case 1
[10954]	118	zmaxu = ze3wu / e3w(ji+1,jj,iku,Kmm)
[5120]	119	! interpolated values of tracers
	120	zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
	121	! gradient of tracers
	122	pgtu(ji,jj,jn) = umask(ji,jj,1) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
	123	ELSE ! case 2
[10954]	124	zmaxu = -ze3wu / e3w(ji,jj,iku,Kmm)
[5120]	125	! interpolated values of tracers
	126	zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
	127	! gradient of tracers
	128	pgtu(ji,jj,jn) = umask(ji,jj,1) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
	129	ENDIF
	130	!
	131	! j- direction
	132	IF( ze3wv >= 0._wp ) THEN ! case 1
[10954]	133	zmaxv = ze3wv / e3w(ji,jj+1,ikv,Kmm)
[5120]	134	! interpolated values of tracers
	135	ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
	136	! gradient of tracers
	137	pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
	138	ELSE ! case 2
[10954]	139	zmaxv = -ze3wv / e3w(ji,jj,ikv,Kmm)
[5120]	140	! interpolated values of tracers
	141	ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
	142	! gradient of tracers
	143	pgtv(ji,jj,jn) = vmask(ji,jj,1) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
	144	ENDIF
	145	END DO
	146	END DO
	147	END DO
[10425]	148	!
	149	CALL lbc_lnk_multi( 'zpshde', pgtu(:,:,:), 'U', -1. , pgtv(:,:,:), 'V', -1. ) ! Lateral boundary cond.
[5836]	150	!
	151	IF( PRESENT( prd ) ) THEN !== horizontal derivative of density anomalies (rd) ==! (optional part)
[7753]	152	pgru(:,:) = 0._wp
	153	pgrv(:,:) = 0._wp ! depth of the partial step level
[5120]	154	DO jj = 1, jpjm1
	155	DO ji = 1, jpim1
	156	iku = mbku(ji,jj)
	157	ikv = mbkv(ji,jj)
[10954]	158	ze3wu = e3w(ji+1,jj ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
	159	ze3wv = e3w(ji ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
	160	IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = gdept(ji ,jj,iku,Kmm) ! i-direction: case 1
	161	ELSE ; zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm) ! - - case 2
[5120]	162	ENDIF
[10954]	163	IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = gdept(ji,jj ,ikv,Kmm) ! j-direction: case 1
	164	ELSE ; zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm) ! - - case 2
[5120]	165	ENDIF
	166	END DO
	167	END DO
[5836]	168	!
	169	CALL eos( zti, zhi, zri ) ! interpolated density from zti, ztj
	170	CALL eos( ztj, zhj, zrj ) ! at the partial step depth output in zri, zrj
	171	!
	172	DO jj = 1, jpjm1 ! Gradient of density at the last level
[5120]	173	DO ji = 1, jpim1
	174	iku = mbku(ji,jj)
	175	ikv = mbkv(ji,jj)
[10954]	176	ze3wu = e3w(ji+1,jj ,iku,Kmm) - e3w(ji,jj,iku,Kmm)
	177	ze3wv = e3w(ji ,jj+1,ikv,Kmm) - e3w(ji,jj,ikv,Kmm)
[5120]	178	IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = umask(ji,jj,1) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1
	179	ELSE ; pgru(ji,jj) = umask(ji,jj,1) * ( prd(ji+1,jj,iku) - zri(ji,jj ) ) ! i: 2
	180	ENDIF
	181	IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = vmask(ji,jj,1) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1
	182	ELSE ; pgrv(ji,jj) = vmask(ji,jj,1) * ( prd(ji,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2
	183	ENDIF
	184	END DO
	185	END DO
[10425]	186	CALL lbc_lnk_multi( 'zpshde', pgru , 'U', -1. , pgrv , 'V', -1. ) ! Lateral boundary conditions
[5120]	187	!
	188	END IF
	189	!
[9019]	190	IF( ln_timing ) CALL timing_stop( 'zps_hde')
[5120]	191	!
	192	END SUBROUTINE zps_hde
[9019]	193
	194
[10954]	195	SUBROUTINE zps_hde_isf( kt, Kmm, kjpt, pta, pgtu, pgtv, pgtui, pgtvi, &
[6140]	196	& prd, pgru, pgrv, pgrui, pgrvi )
[3]	197	!!----------------------------------------------------------------------
[6140]	198	!! * ROUTINE zps_hde_isf *
[3]	199	!!
[2528]	200	!! ** Purpose : Compute the horizontal derivative of T, S and rho
[3]	201	!! at u- and v-points with a linear interpolation for z-coordinate
[6140]	202	!! with partial steps for top (ice shelf) and bottom.
[3]	203	!!
	204	!! ** Method : In z-coord with partial steps, scale factors on last
	205	!! levels are different for each grid point, so that T, S and rd
	206	!! points are not at the same depth as in z-coord. To have horizontal
[6140]	207	!! gradients again, we interpolate T and S at the good depth :
	208	!! For the bottom case:
[3]	209	!! Linear interpolation of T, S
	210	!! Computation of di(tb) and dj(tb) by vertical interpolation:
	211	!! di(t) = t~ - t(i,j,k) or t(i+1,j,k) - t~
	212	!! dj(t) = t~ - t(i,j,k) or t(i,j+1,k) - t~
	213	!! This formulation computes the two cases:
	214	!! CASE 1 CASE 2
	215	!! k-1 ___ ___________ k-1 ___ ___________
	216	!! Ti T~ T~ Ti+1
	217	!! _____ _____
	218	!! k \| \|Ti+1 k Ti \| \|
	219	!! \| \|____ ____\| \|
	220	!! ___ \| \| \| ___ \| \| \|
	221	!!
	222	!! case 1-> e3w(i+1) >= e3w(i) ( and e3w(j+1) >= e3w(j) ) then
	223	!! t~ = t(i+1,j ,k) + (e3w(i+1) - e3w(i)) * dk(Ti+1)/e3w(i+1)
	224	!! ( t~ = t(i ,j+1,k) + (e3w(j+1) - e3w(j)) * dk(Tj+1)/e3w(j+1) )
	225	!! or
	226	!! case 2-> e3w(i+1) <= e3w(i) ( and e3w(j+1) <= e3w(j) ) then
	227	!! t~ = t(i,j,k) + (e3w(i) - e3w(i+1)) * dk(Ti)/e3w(i )
	228	!! ( t~ = t(i,j,k) + (e3w(j) - e3w(j+1)) * dk(Tj)/e3w(j ) )
	229	!! Idem for di(s) and dj(s)
	230	!!
[4990]	231	!! For rho, we call eos which will compute rd~(t~,s~) at the right
	232	!! depth zh from interpolated T and S for the different formulations
	233	!! of the equation of state (eos).
[3]	234	!! Gradient formulation for rho :
[4990]	235	!! di(rho) = rd~ - rd(i,j,k) or rd(i+1,j,k) - rd~
[3]	236	!!
[6140]	237	!! For the top case (ice shelf): As for the bottom case but upside down
	238	!!
[4990]	239	!! ** Action : compute for top and bottom interfaces
[5120]	240	!! - pgtu, pgtv, pgtui, pgtvi: horizontal gradient of tracer at u- & v-points
	241	!! - pgru, pgrv, pgrui, pgtvi: horizontal gradient of rho (if present) at u- & v-points
[2528]	242	!!----------------------------------------------------------------------
[6140]	243	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[10954]	244	INTEGER , INTENT(in ) :: Kmm ! ocean time level index
[6140]	245	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	246	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pta ! 4D tracers fields
	247	REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtu, pgtv ! hor. grad. of ptra at u- & v-pts
	248	REAL(wp), DIMENSION(jpi,jpj, kjpt), INTENT( out) :: pgtui, pgtvi ! hor. grad. of stra at u- & v-pts (ISF)
	249	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ), OPTIONAL :: prd ! 3D density anomaly fields
	250	REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgru, pgrv ! hor. grad of prd at u- & v-pts (bottom)
	251	REAL(wp), DIMENSION(jpi,jpj ), INTENT( out), OPTIONAL :: pgrui, pgrvi ! hor. grad of prd at u- & v-pts (top)
[2715]	252	!
[2528]	253	INTEGER :: ji, jj, jn ! Dummy loop indices
[4990]	254	INTEGER :: iku, ikv, ikum1, ikvm1,ikup1, ikvp1 ! partial step level (ocean bottom level) at u- and v-points
[6140]	255	REAL(wp) :: ze3wu, ze3wv, zmaxu, zmaxv ! temporary scalars
[4990]	256	REAL(wp), DIMENSION(jpi,jpj) :: zri, zrj, zhi, zhj ! NB: 3rd dim=1 to use eos
	257	REAL(wp), DIMENSION(jpi,jpj,kjpt) :: zti, ztj !
[3]	258	!!----------------------------------------------------------------------
[3294]	259	!
[9019]	260	IF( ln_timing ) CALL timing_start( 'zps_hde_isf')
[3294]	261	!
[5836]	262	pgtu (:,:,:) = 0._wp ; pgtv (:,:,:) =0._wp
	263	pgtui(:,:,:) = 0._wp ; pgtvi(:,:,:) =0._wp
	264	zti (:,:,:) = 0._wp ; ztj (:,:,:) =0._wp
	265	zhi (:,: ) = 0._wp ; zhj (:,: ) =0._wp
[3294]	266	!
[2528]	267	DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==!
	268	!
[3]	269	DO jj = 1, jpjm1
[2528]	270	DO ji = 1, jpim1
[6140]	271
	272	iku = mbku(ji,jj); ikum1 = MAX( iku - 1 , 1 ) ! last and before last ocean level at u- & v-points
	273	ikv = mbkv(ji,jj); ikvm1 = MAX( ikv - 1 , 1 ) ! if level first is a p-step, ik.m1=1
[10954]	274	ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
	275	ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
[2528]	276	!
	277	! i- direction
	278	IF( ze3wu >= 0._wp ) THEN ! case 1
[10954]	279	zmaxu = ze3wu / e3w(ji+1,jj,iku,Kmm)
[2528]	280	! interpolated values of tracers
[4990]	281	zti (ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikum1,jn) - pta(ji+1,jj,iku,jn) )
[2528]	282	! gradient of tracers
[6140]	283	pgtu(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
[2528]	284	ELSE ! case 2
[10954]	285	zmaxu = -ze3wu / e3w(ji,jj,iku,Kmm)
[2528]	286	! interpolated values of tracers
[4990]	287	zti (ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikum1,jn) - pta(ji,jj,iku,jn) )
[2528]	288	! gradient of tracers
[6140]	289	pgtu(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
[2528]	290	ENDIF
	291	!
	292	! j- direction
	293	IF( ze3wv >= 0._wp ) THEN ! case 1
[10954]	294	zmaxv = ze3wv / e3w(ji,jj+1,ikv,Kmm)
[2528]	295	! interpolated values of tracers
[4990]	296	ztj (ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvm1,jn) - pta(ji,jj+1,ikv,jn) )
[2528]	297	! gradient of tracers
[6140]	298	pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
[2528]	299	ELSE ! case 2
[10954]	300	zmaxv = -ze3wv / e3w(ji,jj,ikv,Kmm)
[2528]	301	! interpolated values of tracers
[4990]	302	ztj (ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvm1,jn) - pta(ji,jj,ikv,jn) )
[2528]	303	! gradient of tracers
[6140]	304	pgtv(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
[2528]	305	ENDIF
[6140]	306
[3]	307	END DO
	308	END DO
[2528]	309	END DO
[10425]	310	!
	311	CALL lbc_lnk_multi( 'zpshde', pgtu(:,:,:), 'U', -1. , pgtv(:,:,:), 'V', -1. ) ! Lateral boundary cond.
[3]	312
[6140]	313	! horizontal derivative of density anomalies (rd)
	314	IF( PRESENT( prd ) ) THEN ! depth of the partial step level
	315	pgru(:,:)=0.0_wp ; pgrv(:,:)=0.0_wp ;
[5836]	316	!
[6140]	317	DO jj = 1, jpjm1
[2528]	318	DO ji = 1, jpim1
[6140]	319
[2528]	320	iku = mbku(ji,jj)
	321	ikv = mbkv(ji,jj)
[10954]	322	ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
	323	ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
[5836]	324	!
[10954]	325	IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = gdept(ji ,jj,iku,Kmm) ! i-direction: case 1
	326	ELSE ; zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm) ! - - case 2
[2528]	327	ENDIF
[10954]	328	IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = gdept(ji,jj ,ikv,Kmm) ! j-direction: case 1
	329	ELSE ; zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm) ! - - case 2
[2528]	330	ENDIF
[6140]	331
[2528]	332	END DO
[3]	333	END DO
	334
[6140]	335	! Compute interpolated rd from zti, ztj for the 2 cases at the depth of the partial
	336	! step and store it in zri, zrj for each case
	337	CALL eos( zti, zhi, zri )
	338	CALL eos( ztj, zhj, zrj )
	339
[5836]	340	DO jj = 1, jpjm1 ! Gradient of density at the last level
[4990]	341	DO ji = 1, jpim1
[6140]	342	iku = mbku(ji,jj)
	343	ikv = mbkv(ji,jj)
[10954]	344	ze3wu = gdept(ji+1,jj,iku,Kmm) - gdept(ji,jj,iku,Kmm)
	345	ze3wv = gdept(ji,jj+1,ikv,Kmm) - gdept(ji,jj,ikv,Kmm)
[6140]	346
	347	IF( ze3wu >= 0._wp ) THEN ; pgru(ji,jj) = ssumask(ji,jj) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1
	348	ELSE ; pgru(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj,iku) - zri(ji,jj ) ) ! i: 2
[4990]	349	ENDIF
[6140]	350	IF( ze3wv >= 0._wp ) THEN ; pgrv(ji,jj) = ssvmask(ji,jj) * ( zrj(ji,jj ) - prd(ji,jj,ikv) ) ! j: 1
	351	ELSE ; pgrv(ji,jj) = ssvmask(ji,jj) * ( prd(ji,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2
[4990]	352	ENDIF
[6140]	353
[4990]	354	END DO
	355	END DO
[6140]	356
[10425]	357	CALL lbc_lnk_multi( 'zpshde', pgru , 'U', -1. , pgrv , 'V', -1. ) ! Lateral boundary conditions
[4990]	358	!
	359	END IF
[5836]	360	!
	361	! !== (ISH) compute grui and gruvi ==!
	362	!
[4990]	363	DO jn = 1, kjpt !== Interpolation of tracers at the last ocean level ==! !
	364	DO jj = 1, jpjm1
	365	DO ji = 1, jpim1
[6140]	366	iku = miku(ji,jj); ikup1 = miku(ji,jj) + 1
	367	ikv = mikv(ji,jj); ikvp1 = mikv(ji,jj) + 1
[4990]	368	!
	369	! (ISF) case partial step top and bottom in adjacent cell in vertical
	370	! cannot used e3w because if 2 cell water column, we have ps at top and bottom
	371	! in this case e3w(i,j) - e3w(i,j+1) is not the distance between Tj~ and Tj
	372	! the only common depth between cells (i,j) and (i,j+1) is gdepw_0
[10954]	373	ze3wu = gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
	374	ze3wv = gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm)
[6140]	375
[4990]	376	! i- direction
	377	IF( ze3wu >= 0._wp ) THEN ! case 1
[10954]	378	zmaxu = ze3wu / e3w(ji+1,jj,ikup1,Kmm)
[4990]	379	! interpolated values of tracers
[6140]	380	zti(ji,jj,jn) = pta(ji+1,jj,iku,jn) + zmaxu * ( pta(ji+1,jj,ikup1,jn) - pta(ji+1,jj,iku,jn) )
[4990]	381	! gradient of tracers
[6140]	382	pgtui(ji,jj,jn) = ssumask(ji,jj) * ( zti(ji,jj,jn) - pta(ji,jj,iku,jn) )
[4990]	383	ELSE ! case 2
[10954]	384	zmaxu = - ze3wu / e3w(ji,jj,ikup1,Kmm)
[4990]	385	! interpolated values of tracers
[6140]	386	zti(ji,jj,jn) = pta(ji,jj,iku,jn) + zmaxu * ( pta(ji,jj,ikup1,jn) - pta(ji,jj,iku,jn) )
[4990]	387	! gradient of tracers
[6140]	388	pgtui(ji,jj,jn) = ssumask(ji,jj) * ( pta(ji+1,jj,iku,jn) - zti(ji,jj,jn) )
[4990]	389	ENDIF
	390	!
	391	! j- direction
	392	IF( ze3wv >= 0._wp ) THEN ! case 1
[10954]	393	zmaxv = ze3wv / e3w(ji,jj+1,ikvp1,Kmm)
[4990]	394	! interpolated values of tracers
[6140]	395	ztj(ji,jj,jn) = pta(ji,jj+1,ikv,jn) + zmaxv * ( pta(ji,jj+1,ikvp1,jn) - pta(ji,jj+1,ikv,jn) )
[4990]	396	! gradient of tracers
[6140]	397	pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( ztj(ji,jj,jn) - pta(ji,jj,ikv,jn) )
[4990]	398	ELSE ! case 2
[10954]	399	zmaxv = - ze3wv / e3w(ji,jj,ikvp1,Kmm)
[4990]	400	! interpolated values of tracers
[6140]	401	ztj(ji,jj,jn) = pta(ji,jj,ikv,jn) + zmaxv * ( pta(ji,jj,ikvp1,jn) - pta(ji,jj,ikv,jn) )
[4990]	402	! gradient of tracers
[6140]	403	pgtvi(ji,jj,jn) = ssvmask(ji,jj) * ( pta(ji,jj+1,ikv,jn) - ztj(ji,jj,jn) )
[4990]	404	ENDIF
[6140]	405
	406	END DO
	407	END DO
[4990]	408	!
	409	END DO
[10425]	410	CALL lbc_lnk_multi( 'zpshde', pgtui(:,:,:), 'U', -1. , pgtvi(:,:,:), 'V', -1. ) ! Lateral boundary cond.
[4990]	411
[5836]	412	IF( PRESENT( prd ) ) THEN !== horizontal derivative of density anomalies (rd) ==! (optional part)
	413	!
[6140]	414	pgrui(:,:) =0.0_wp; pgrvi(:,:) =0.0_wp;
	415	DO jj = 1, jpjm1
[4990]	416	DO ji = 1, jpim1
[6140]	417
[4990]	418	iku = miku(ji,jj)
	419	ikv = mikv(ji,jj)
[10954]	420	ze3wu = gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
	421	ze3wv = gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm)
[5836]	422	!
[10954]	423	IF( ze3wu >= 0._wp ) THEN ; zhi(ji,jj) = gdept(ji ,jj,iku,Kmm) ! i-direction: case 1
	424	ELSE ; zhi(ji,jj) = gdept(ji+1,jj,iku,Kmm) ! - - case 2
[4990]	425	ENDIF
[6140]	426
[10954]	427	IF( ze3wv >= 0._wp ) THEN ; zhj(ji,jj) = gdept(ji,jj ,ikv,Kmm) ! j-direction: case 1
	428	ELSE ; zhj(ji,jj) = gdept(ji,jj+1,ikv,Kmm) ! - - case 2
[4990]	429	ENDIF
[6140]	430
[4990]	431	END DO
	432	END DO
[5836]	433	!
	434	CALL eos( zti, zhi, zri ) ! interpolated density from zti, ztj
	435	CALL eos( ztj, zhj, zrj ) ! at the partial step depth output in zri, zrj
	436	!
	437	DO jj = 1, jpjm1 ! Gradient of density at the last level
[2528]	438	DO ji = 1, jpim1
[6140]	439	iku = miku(ji,jj)
	440	ikv = mikv(ji,jj)
[10954]	441	ze3wu = gdept(ji,jj,iku,Kmm) - gdept(ji+1,jj,iku,Kmm)
	442	ze3wv = gdept(ji,jj,ikv,Kmm) - gdept(ji,jj+1,ikv,Kmm)
[6140]	443
	444	IF( ze3wu >= 0._wp ) THEN ; pgrui(ji,jj) = ssumask(ji,jj) * ( zri(ji ,jj ) - prd(ji,jj,iku) ) ! i: 1
	445	ELSE ; pgrui(ji,jj) = ssumask(ji,jj) * ( prd(ji+1,jj ,iku) - zri(ji,jj ) ) ! i: 2
[2528]	446	ENDIF
[6140]	447	IF( ze3wv >= 0._wp ) THEN ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( zrj(ji ,jj ) - prd(ji,jj,ikv) ) ! j: 1
	448	ELSE ; pgrvi(ji,jj) = ssvmask(ji,jj) * ( prd(ji ,jj+1,ikv) - zrj(ji,jj ) ) ! j: 2
[2528]	449	ENDIF
[6140]	450
[2528]	451	END DO
[3]	452	END DO
[10425]	453	CALL lbc_lnk_multi( 'zpshde', pgrui, 'U', -1. , pgrvi, 'V', -1. ) ! Lateral boundary conditions
[2528]	454	!
[4990]	455	END IF
[2528]	456	!
[9019]	457	IF( ln_timing ) CALL timing_stop( 'zps_hde_isf')
[2715]	458	!
[5120]	459	END SUBROUTINE zps_hde_isf
[9019]	460
[3]	461	!!======================================================================
	462	END MODULE zpshde

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

Download in other formats: