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

Context Navigation

source: NEMO/trunk/src/OCE/FLO/flo4rk.F90 @ 12489

Last change on this file since 12489 was 12489, checked in by davestorkey, 4 years ago

Preparation for new timestepping scheme #2390.
Main changes:

Initial euler timestep now handled in stp and not in TRA/DYN routines.
Renaming of all timestep parameters. In summary, the namelist parameter is now rn_Dt and the current timestep is rDt (and rDt_ice, rDt_trc etc).
Renaming of a few miscellaneous parameters, eg. atfp -> rn_atfp (namelist parameter used everywhere) and rau0 -> rho0.

This version gives bit-comparable results to the previous version of the trunk.

Property svn:keywords set to Id

File size: 18.5 KB

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

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

Download in other formats: