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/trunk/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/TRA/zpshde.F90 @ 14719

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

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

Download in other formats: