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 trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/zpshde.F90 @ 7702

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

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

Download in other formats: