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/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13723_KERNEL-01_Amy_Mike_newHPGschemes/src/OCE/TRA/zpshde.F90 @ 14046

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

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

Download in other formats: