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.

flo4rk.F90 in branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/FLO – NEMO

Context Navigation

source: branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/FLO/flo4rk.F90 @ 5591

Last change on this file since 5591 was 5477, checked in by cguiavarch, 9 years ago
Clear svn keywords from UKMO/dev_r5107_hadgem3_cplseq
File size: 19.8 KB

Rev	Line
[3]	1	MODULE flo4rk
	2	!!======================================================================
	3	!! * MODULE flo4rk *
	4	!! Ocean floats : trajectory computation using a 4th order Runge-Kutta
	5	!!======================================================================
	6	#if defined key_floats \|\| defined key_esopa
	7	!!----------------------------------------------------------------------
	8	!! 'key_floats' float trajectories
	9	!!----------------------------------------------------------------------
	10	!! flo_4rk : Compute the geographical position of floats
	11	!! flo_interp : interpolation
	12	!!----------------------------------------------------------------------
	13	USE flo_oce ! ocean drifting floats
	14	USE oce ! ocean dynamics and tracers
	15	USE dom_oce ! ocean space and time domain
[16]	16	USE in_out_manager ! I/O manager
[3294]	17	USE wrk_nemo ! working array
[3]	18
	19	IMPLICIT NONE
	20	PRIVATE
	21
[2528]	22	PUBLIC flo_4rk ! routine called by floats.F90
[3]	23
[2528]	24	! ! RK4 and Lagrange interpolation coefficients
	25	REAL(wp), DIMENSION (4) :: tcoef1 = (/ 1.0 , 0.5 , 0.5 , 0.0 /) !
	26	REAL(wp), DIMENSION (4) :: tcoef2 = (/ 0.0 , 0.5 , 0.5 , 1.0 /) !
	27	REAL(wp), DIMENSION (4) :: scoef2 = (/ 1.0 , 2.0 , 2.0 , 1.0 /) !
	28	REAL(wp), DIMENSION (4) :: rcoef = (/-1./6. , 1./2. ,-1./2. , 1./6. /) !
	29	REAL(wp), DIMENSION (3) :: scoef1 = (/ 0.5 , 0.5 , 1.0 /) !
	30
	31	!! * Substitutions
	32	# include "domzgr_substitute.h90"
[3]	33	!!----------------------------------------------------------------------
[2528]	34	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[5477]	35	!! $Id$
[2528]	36	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	37	!!----------------------------------------------------------------------
	38	CONTAINS
	39
	40	SUBROUTINE flo_4rk( kt )
	41	!!----------------------------------------------------------------------
	42	!! * ROUTINE flo_4rk *
	43	!!
	44	!! ** Purpose : Compute the geographical position (lat,lon,depth)
	45	!! of each float at each time step.
	46	!!
	47	!! ** Method : The position of a float is computed with a 4th order
	48	!! Runge-Kutta scheme and and Lagrange interpolation.
	49	!! We need to know the velocity field, the old positions of the
	50	!! floats and the grid defined on the domain.
[2528]	51	!!----------------------------------------------------------------------
	52	INTEGER, INTENT(in) :: kt ! ocean time-step index
[3]	53	!!
	54	INTEGER :: jfl, jind ! dummy loop indices
[3294]	55	INTEGER :: ierror ! error value
	56
	57	REAL(wp), POINTER, DIMENSION(:) :: zgifl , zgjfl , zgkfl ! index RK positions
	58	REAL(wp), POINTER, DIMENSION(:) :: zufl , zvfl , zwfl ! interpolated velocity at the float position
	59	REAL(wp), POINTER, DIMENSION(:,:) :: zrkxfl, zrkyfl, zrkzfl ! RK coefficients
[3]	60	!!---------------------------------------------------------------------
[3294]	61	CALL wrk_alloc( jpnfl, zgifl , zgjfl , zgkfl , zufl, zvfl, zwfl)
	62	CALL wrk_alloc( jpnfl, 4, zrkxfl, zrkyfl, zrkzfl )
	63	!
	64	IF( ierror /= 0 ) THEN
	65	WRITE(numout,*) 'flo_4rk: allocation of workspace arrays failed'
	66	ENDIF
	67
[3]	68
	69	IF( kt == nit000 ) THEN
	70	IF(lwp) WRITE(numout,*)
	71	IF(lwp) WRITE(numout,*) 'flo_4rk : compute Runge Kutta trajectories for floats '
	72	IF(lwp) WRITE(numout,*) '~~~~~~~'
	73	ENDIF
	74
	75	! Verification of the floats positions. If one of them leave the domain
	76	! domain we replace the float near the border.
	77	DO jfl = 1, jpnfl
	78	! i-direction
	79	IF( tpifl(jfl) <= 1.5 ) THEN
	80	IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
	81	IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the WEST border.'
	82	tpifl(jfl) = tpifl(jfl) + 1.
	83	IF(lwp)WRITE(numout,*)'New initialisation for this float at i=',tpifl(jfl)
	84	ENDIF
	85
	86	IF( tpifl(jfl) >= jpi-.5 ) THEN
	87	IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
	88	IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the EAST border.'
	89	tpifl(jfl) = tpifl(jfl) - 1.
	90	IF(lwp)WRITE(numout,*)'New initialisation for this float at i=', tpifl(jfl)
	91	ENDIF
	92	! j-direction
	93	IF( tpjfl(jfl) <= 1.5 ) THEN
	94	IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
	95	IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the SOUTH border.'
	96	tpjfl(jfl) = tpjfl(jfl) + 1.
	97	IF(lwp)WRITE(numout,*)'New initialisation for this float at j=', tpjfl(jfl)
	98	ENDIF
	99
	100	IF( tpjfl(jfl) >= jpj-.5 ) THEN
	101	IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
	102	IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the NORTH border.'
	103	tpjfl(jfl) = tpjfl(jfl) - 1.
	104	IF(lwp)WRITE(numout,*)'New initialisation for this float at j=', tpjfl(jfl)
	105	ENDIF
	106	! k-direction
	107	IF( tpkfl(jfl) <= .5 ) THEN
	108	IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
	109	IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the TOP border.'
	110	tpkfl(jfl) = tpkfl(jfl) + 1.
	111	IF(lwp)WRITE(numout,*)'New initialisation for this float at k=', tpkfl(jfl)
	112	ENDIF
	113
	114	IF( tpkfl(jfl) >= jpk-.5 ) THEN
	115	IF(lwp)WRITE(numout,*)'!!!!!!!!!!!!! WARNING !!!!!!!!!!!!!!!!'
	116	IF(lwp)WRITE(numout,*)'The float',jfl,'is out of the domain at the BOTTOM border.'
	117	tpkfl(jfl) = tpkfl(jfl) - 1.
	118	IF(lwp)WRITE(numout,*)'New initialisation for this float at k=', tpkfl(jfl)
	119	ENDIF
	120	END DO
	121
	122	! 4 steps of Runge-Kutta algorithme
	123	! initialisation of the positions
	124
	125	DO jfl = 1, jpnfl
	126	zgifl(jfl) = tpifl(jfl)
	127	zgjfl(jfl) = tpjfl(jfl)
	128	zgkfl(jfl) = tpkfl(jfl)
	129	END DO
	130
[2528]	131	DO jind = 1, 4
	132
[3]	133	! for each step we compute the compute the velocity with Lagrange interpolation
[2528]	134	CALL flo_interp( zgifl, zgjfl, zgkfl, zufl, zvfl, zwfl, jind )
[3]	135
	136	! computation of Runge-Kutta factor
	137	DO jfl = 1, jpnfl
	138	zrkxfl(jfl,jind) = rdt*zufl(jfl)
	139	zrkyfl(jfl,jind) = rdt*zvfl(jfl)
	140	zrkzfl(jfl,jind) = rdt*zwfl(jfl)
	141	END DO
	142	IF( jind /= 4 ) THEN
	143	DO jfl = 1, jpnfl
	144	zgifl(jfl) = (tpifl(jfl)) + scoef1(jind)*zrkxfl(jfl,jind)
	145	zgjfl(jfl) = (tpjfl(jfl)) + scoef1(jind)*zrkyfl(jfl,jind)
	146	zgkfl(jfl) = (tpkfl(jfl)) + scoef1(jind)*zrkzfl(jfl,jind)
	147	END DO
	148	ENDIF
	149	END DO
	150	DO jind = 1, 4
	151	DO jfl = 1, jpnfl
	152	tpifl(jfl) = tpifl(jfl) + scoef2(jind)*zrkxfl(jfl,jind)/6.
	153	tpjfl(jfl) = tpjfl(jfl) + scoef2(jind)*zrkyfl(jfl,jind)/6.
	154	tpkfl(jfl) = tpkfl(jfl) + scoef2(jind)*zrkzfl(jfl,jind)/6.
	155	END DO
	156	END DO
[2528]	157	!
[3294]	158	CALL wrk_dealloc( jpnfl, zgifl , zgjfl , zgkfl , zufl, zvfl, zwfl)
	159	CALL wrk_dealloc( jpnfl, 4, zrkxfl, zrkyfl, zrkzfl )
	160	!
[3]	161	END SUBROUTINE flo_4rk
	162
	163
	164	SUBROUTINE flo_interp( pxt , pyt , pzt , &
[2528]	165	& pufl, pvfl, pwfl, ki )
[3]	166	!!----------------------------------------------------------------------
	167	!! * ROUTINE flointerp *
	168	!!
	169	!! ** Purpose : Interpolation of the velocity on the float position
	170	!!
	171	!! ** Method : Lagrange interpolation with the 64 neighboring
	172	!! points. This routine is call 4 time at each time step to
	173	!! compute velocity at the date and the position we need to
	174	!! integrated with RK method.
[2528]	175	!!----------------------------------------------------------------------
	176	REAL(wp) , DIMENSION(jpnfl), INTENT(in ) :: pxt , pyt , pzt ! position of the float
	177	REAL(wp) , DIMENSION(jpnfl), INTENT( out) :: pufl, pvfl, pwfl ! velocity at this position
	178	INTEGER , INTENT(in ) :: ki !
[3]	179	!!
[2528]	180	INTEGER :: jfl, jind1, jind2, jind3 ! dummy loop indices
	181	REAL(wp) :: zsumu, zsumv, zsumw ! local scalar
[3294]	182	INTEGER , POINTER, DIMENSION(:) :: iilu, ijlu, iklu ! nearest neighbour INDEX-u
	183	INTEGER , POINTER, DIMENSION(:) :: iilv, ijlv, iklv ! nearest neighbour INDEX-v
	184	INTEGER , POINTER, DIMENSION(:) :: iilw, ijlw, iklw ! nearest neighbour INDEX-w
	185	INTEGER , POINTER, DIMENSION(:,:) :: iidu, ijdu, ikdu ! 64 nearest neighbour INDEX-u
	186	INTEGER , POINTER, DIMENSION(:,:) :: iidv, ijdv, ikdv ! 64 nearest neighbour INDEX-v
	187	INTEGER , POINTER, DIMENSION(:,:) :: iidw, ijdw, ikdw ! 64 nearest neighbour INDEX-w
	188	REAL(wp) , POINTER, DIMENSION(:,:) :: zlagxu, zlagyu, zlagzu ! Lagrange coefficients
	189	REAL(wp) , POINTER, DIMENSION(:,:) :: zlagxv, zlagyv, zlagzv ! - -
	190	REAL(wp) , POINTER, DIMENSION(:,:) :: zlagxw, zlagyw, zlagzw ! - -
	191	REAL(wp) , POINTER, DIMENSION(:,:,:,:) :: ztufl , ztvfl , ztwfl ! velocity at choosen time step
[3]	192	!!---------------------------------------------------------------------
[3294]	193	CALL wrk_alloc( jpnfl, iilu, ijlu, iklu, iilv, ijlv, iklv, iilw, ijlw, iklw )
	194	CALL wrk_alloc( jpnfl, 4, iidu, ijdu, ikdu, iidv, ijdv, ikdv, iidw, ijdw, ikdw )
	195	CALL wrk_alloc( jpnfl, 4, zlagxu, zlagyu, zlagzu, zlagxv, zlagyv, zlagzv, zlagxw, zlagyw, zlagzw )
	196	CALL wrk_alloc( jpnfl, 4, 4, 4, ztufl , ztvfl , ztwfl )
	197
[3]	198	! Interpolation of U velocity
	199
	200	! nearest neightboring point for computation of u
	201	DO jfl = 1, jpnfl
	202	iilu(jfl) = INT(pxt(jfl)-.5)
	203	ijlu(jfl) = INT(pyt(jfl)-.5)
	204	iklu(jfl) = INT(pzt(jfl))
	205	END DO
	206
	207	! 64 neightboring points for computation of u
	208	DO jind1 = 1, 4
	209	DO jfl = 1, jpnfl
	210	! i-direction
[2528]	211	IF( iilu(jfl) <= 2 ) THEN ; iidu(jfl,jind1) = jind1
[3]	212	ELSE
[2528]	213	IF( iilu(jfl) >= jpi-1 ) THEN ; iidu(jfl,jind1) = jpi + jind1 - 4
	214	ELSE ; iidu(jfl,jind1) = iilu(jfl) + jind1 - 2
[3]	215	ENDIF
	216	ENDIF
	217	! j-direction
[2528]	218	IF( ijlu(jfl) <= 2 ) THEN ; ijdu(jfl,jind1) = jind1
[3]	219	ELSE
[2528]	220	IF( ijlu(jfl) >= jpj-1 ) THEN ; ijdu(jfl,jind1) = jpj + jind1 - 4
	221	ELSE ; ijdu(jfl,jind1) = ijlu(jfl) + jind1 - 2
[3]	222	ENDIF
	223	ENDIF
	224	! k-direction
[2528]	225	IF( iklu(jfl) <= 2 ) THEN ; ikdu(jfl,jind1) = jind1
[3]	226	ELSE
[2528]	227	IF( iklu(jfl) >= jpk-1 ) THEN ; ikdu(jfl,jind1) = jpk + jind1 - 4
	228	ELSE ; ikdu(jfl,jind1) = iklu(jfl) + jind1 - 2
[3]	229	ENDIF
	230	ENDIF
	231	END DO
	232	END DO
	233
	234	! Lagrange coefficients
	235	DO jfl = 1, jpnfl
	236	DO jind1 = 1, 4
	237	zlagxu(jfl,jind1) = 1.
	238	zlagyu(jfl,jind1) = 1.
	239	zlagzu(jfl,jind1) = 1.
	240	END DO
	241	END DO
	242	DO jind1 = 1, 4
	243	DO jind2 = 1, 4
	244	DO jfl= 1, jpnfl
	245	IF( jind1 /= jind2 ) THEN
	246	zlagxu(jfl,jind1) = zlagxu(jfl,jind1) * ( pxt(jfl)-(float(iidu(jfl,jind2))+.5) )
	247	zlagyu(jfl,jind1) = zlagyu(jfl,jind1) * ( pyt(jfl)-(float(ijdu(jfl,jind2))) )
	248	zlagzu(jfl,jind1) = zlagzu(jfl,jind1) * ( pzt(jfl)-(float(ikdu(jfl,jind2))) )
	249	ENDIF
	250	END DO
	251	END DO
	252	END DO
	253
	254	! velocity when we compute at middle time step
	255
	256	DO jfl = 1, jpnfl
	257	DO jind1 = 1, 4
	258	DO jind2 = 1, 4
	259	DO jind3 = 1, 4
	260	ztufl(jfl,jind1,jind2,jind3) = &
[2528]	261	& ( tcoef1(ki) * ub(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) + &
	262	& tcoef2(ki) * un(iidu(jfl,jind1),ijdu(jfl,jind2),ikdu(jfl,jind3)) ) &
[3]	263	& / e1u(iidu(jfl,jind1),ijdu(jfl,jind2))
	264	END DO
	265	END DO
	266	END DO
	267
	268	zsumu = 0.
	269	DO jind1 = 1, 4
	270	DO jind2 = 1, 4
	271	DO jind3 = 1, 4
	272	zsumu = zsumu + ztufl(jfl,jind1,jind2,jind3) * zlagxu(jfl,jind1) * zlagyu(jfl,jind2) &
	273	& * zlagzu(jfl,jind3) * rcoef(jind1)rcoef(jind2)rcoef(jind3)
	274	END DO
	275	END DO
	276	END DO
	277	pufl(jfl) = zsumu
	278	END DO
	279
	280	! Interpolation of V velocity
	281
	282	! nearest neightboring point for computation of v
	283	DO jfl = 1, jpnfl
	284	iilv(jfl) = INT(pxt(jfl)-.5)
	285	ijlv(jfl) = INT(pyt(jfl)-.5)
	286	iklv(jfl) = INT(pzt(jfl))
	287	END DO
	288
	289	! 64 neightboring points for computation of v
	290	DO jind1 = 1, 4
	291	DO jfl = 1, jpnfl
	292	! i-direction
[2528]	293	IF( iilv(jfl) <= 2 ) THEN ; iidv(jfl,jind1) = jind1
[3]	294	ELSE
[2528]	295	IF( iilv(jfl) >= jpi-1 ) THEN ; iidv(jfl,jind1) = jpi + jind1 - 4
	296	ELSE ; iidv(jfl,jind1) = iilv(jfl) + jind1 - 2
[3]	297	ENDIF
	298	ENDIF
	299	! j-direction
[2528]	300	IF( ijlv(jfl) <= 2 ) THEN ; ijdv(jfl,jind1) = jind1
[3]	301	ELSE
[2528]	302	IF( ijlv(jfl) >= jpj-1 ) THEN ; ijdv(jfl,jind1) = jpj + jind1 - 4
	303	ELSE ; ijdv(jfl,jind1) = ijlv(jfl) + jind1 - 2
[3]	304	ENDIF
	305	ENDIF
	306	! k-direction
[2528]	307	IF( iklv(jfl) <= 2 ) THEN ; ikdv(jfl,jind1) = jind1
[3]	308	ELSE
[2528]	309	IF( iklv(jfl) >= jpk-1 ) THEN ; ikdv(jfl,jind1) = jpk + jind1 - 4
	310	ELSE ; ikdv(jfl,jind1) = iklv(jfl) + jind1 - 2
[3]	311	ENDIF
	312	ENDIF
	313	END DO
	314	END DO
	315
	316	! Lagrange coefficients
	317
	318	DO jfl = 1, jpnfl
	319	DO jind1 = 1, 4
	320	zlagxv(jfl,jind1) = 1.
	321	zlagyv(jfl,jind1) = 1.
	322	zlagzv(jfl,jind1) = 1.
	323	END DO
	324	END DO
	325
	326	DO jind1 = 1, 4
	327	DO jind2 = 1, 4
	328	DO jfl = 1, jpnfl
	329	IF( jind1 /= jind2 ) THEN
[2528]	330	zlagxv(jfl,jind1)= zlagxv(jfl,jind1)*(pxt(jfl) - (float(iidv(jfl,jind2)) ) )
[3]	331	zlagyv(jfl,jind1)= zlagyv(jfl,jind1)*(pyt(jfl) - (float(ijdv(jfl,jind2))+.5) )
[2528]	332	zlagzv(jfl,jind1)= zlagzv(jfl,jind1)*(pzt(jfl) - (float(ikdv(jfl,jind2)) ) )
[3]	333	ENDIF
	334	END DO
	335	END DO
	336	END DO
	337
	338	! velocity when we compute at middle time step
	339
	340	DO jfl = 1, jpnfl
	341	DO jind1 = 1, 4
	342	DO jind2 = 1, 4
	343	DO jind3 = 1 ,4
	344	ztvfl(jfl,jind1,jind2,jind3)= &
[2528]	345	& ( tcoef1(ki) * vb(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3)) + &
	346	& tcoef2(ki) * vn(iidv(jfl,jind1),ijdv(jfl,jind2),ikdv(jfl,jind3)) ) &
[3]	347	& / e2v(iidv(jfl,jind1),ijdv(jfl,jind2))
	348	END DO
	349	END DO
	350	END DO
	351
	352	zsumv=0.
	353	DO jind1 = 1, 4
	354	DO jind2 = 1, 4
	355	DO jind3 = 1, 4
	356	zsumv = zsumv + ztvfl(jfl,jind1,jind2,jind3) * zlagxv(jfl,jind1) * zlagyv(jfl,jind2) &
	357	& * zlagzv(jfl,jind3) * rcoef(jind1)rcoef(jind2)rcoef(jind3)
	358	END DO
	359	END DO
	360	END DO
	361	pvfl(jfl) = zsumv
	362	END DO
	363
	364	! Interpolation of W velocity
	365
	366	! nearest neightboring point for computation of w
	367	DO jfl = 1, jpnfl
[2528]	368	iilw(jfl) = INT( pxt(jfl) )
	369	ijlw(jfl) = INT( pyt(jfl) )
	370	iklw(jfl) = INT( pzt(jfl)+.5)
[3]	371	END DO
	372
	373	! 64 neightboring points for computation of w
	374	DO jind1 = 1, 4
	375	DO jfl = 1, jpnfl
	376	! i-direction
[2528]	377	IF( iilw(jfl) <= 2 ) THEN ; iidw(jfl,jind1) = jind1
[3]	378	ELSE
[2528]	379	IF( iilw(jfl) >= jpi-1 ) THEN ; iidw(jfl,jind1) = jpi + jind1 - 4
	380	ELSE ; iidw(jfl,jind1) = iilw(jfl) + jind1 - 2
[3]	381	ENDIF
	382	ENDIF
	383	! j-direction
[2528]	384	IF( ijlw(jfl) <= 2 ) THEN ; ijdw(jfl,jind1) = jind1
[3]	385	ELSE
[2528]	386	IF( ijlw(jfl) >= jpj-1 ) THEN ; ijdw(jfl,jind1) = jpj + jind1 - 4
	387	ELSE ; ijdw(jfl,jind1) = ijlw(jfl) + jind1 - 2
[3]	388	ENDIF
	389	ENDIF
	390	! k-direction
[2528]	391	IF( iklw(jfl) <= 2 ) THEN ; ikdw(jfl,jind1) = jind1
[3]	392	ELSE
[2528]	393	IF( iklw(jfl) >= jpk-1 ) THEN ; ikdw(jfl,jind1) = jpk + jind1 - 4
	394	ELSE ; ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
[3]	395	ENDIF
	396	ENDIF
	397	END DO
	398	END DO
	399	DO jind1 = 1, 4
	400	DO jfl = 1, jpnfl
[2528]	401	IF( iklw(jfl) <= 2 ) THEN ; ikdw(jfl,jind1) = jind1
[3]	402	ELSE
[2528]	403	IF( iklw(jfl) >= jpk-1 ) THEN ; ikdw(jfl,jind1) = jpk + jind1 - 4
	404	ELSE ; ikdw(jfl,jind1) = iklw(jfl) + jind1 - 2
[3]	405	ENDIF
	406	ENDIF
	407	END DO
	408	END DO
	409
	410	! Lagrange coefficients for w interpolation
	411	DO jfl = 1, jpnfl
	412	DO jind1 = 1, 4
	413	zlagxw(jfl,jind1) = 1.
	414	zlagyw(jfl,jind1) = 1.
	415	zlagzw(jfl,jind1) = 1.
	416	END DO
	417	END DO
	418	DO jind1 = 1, 4
	419	DO jind2 = 1, 4
	420	DO jfl = 1, jpnfl
	421	IF( jind1 /= jind2 ) THEN
[2528]	422	zlagxw(jfl,jind1) = zlagxw(jfl,jind1) * (pxt(jfl) - (float(iidw(jfl,jind2)) ) )
	423	zlagyw(jfl,jind1) = zlagyw(jfl,jind1) * (pyt(jfl) - (float(ijdw(jfl,jind2)) ) )
[3]	424	zlagzw(jfl,jind1) = zlagzw(jfl,jind1) * (pzt(jfl) - (float(ikdw(jfl,jind2))-.5) )
	425	ENDIF
	426	END DO
	427	END DO
	428	END DO
	429
	430	! velocity w when we compute at middle time step
	431	DO jfl = 1, jpnfl
	432	DO jind1 = 1, 4
	433	DO jind2 = 1, 4
	434	DO jind3 = 1, 4
	435	ztwfl(jfl,jind1,jind2,jind3)= &
[2528]	436	& ( tcoef1(ki) * wb(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3))+ &
	437	& tcoef2(ki) * wn(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3)) ) &
	438	& / fse3w(iidw(jfl,jind1),ijdw(jfl,jind2),ikdw(jfl,jind3))
[3]	439	END DO
	440	END DO
	441	END DO
	442
[2528]	443	zsumw = 0.e0
[3]	444	DO jind1 = 1, 4
	445	DO jind2 = 1, 4
	446	DO jind3 = 1, 4
	447	zsumw = zsumw + ztwfl(jfl,jind1,jind2,jind3) * zlagxw(jfl,jind1) * zlagyw(jfl,jind2) &
	448	& * zlagzw(jfl,jind3) * rcoef(jind1)rcoef(jind2)rcoef(jind3)
	449	END DO
	450	END DO
	451	END DO
	452	pwfl(jfl) = zsumw
	453	END DO
[2528]	454	!
[3294]	455	CALL wrk_dealloc( jpnfl, iilu, ijlu, iklu, iilv, ijlv, iklv, iilw, ijlw, iklw )
	456	CALL wrk_dealloc( jpnfl, 4, iidu, ijdu, ikdu, iidv, ijdv, ikdv, iidw, ijdw, ikdw )
	457	CALL wrk_dealloc( jpnfl, 4, zlagxu, zlagyu, zlagzu, zlagxv, zlagyv, zlagzv, zlagxw, zlagyw, zlagzw )
	458	CALL wrk_dealloc( jpnfl, 4, 4, 4, ztufl , ztvfl , ztwfl )
	459	!
[3]	460	END SUBROUTINE flo_interp
	461
	462	# else
	463	!!----------------------------------------------------------------------
[2528]	464	!! No floats Dummy module
[3]	465	!!----------------------------------------------------------------------
	466	#endif
	467
	468	!!======================================================================
	469	END MODULE flo4rk

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

Download in other formats: