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.

stprk3.F90 in NEMO/trunk/src/SWE – NEMO

Context Navigation

source: NEMO/trunk/src/SWE/stprk3.F90 @ 14239

Last change on this file since 14239 was 14239, checked in by smasson, 3 years ago
trunk: replace key_iomput by key_xios
File size: 16.7 KB

Rev	Line
[13769]	1	MODULE stprk3
[12983]	2	!!======================================================================
[13769]	3	!! * MODULE stprk3 *
[12983]	4	!! Time-stepping : manager of the shallow water equation time stepping
[13604]	5	!! 3rd order Runge-Kutta scheme
[12983]	6	!!======================================================================
	7	!! History : NEMO ! 2020-03 (A. Nasser, G. Madec) Original code from 4.0.2
[13604]	8	!! - ! 2020-10 (S. Techene, G. Madec) cleanning
[12983]	9	!!----------------------------------------------------------------------
	10
	11	!!----------------------------------------------------------------------
[13604]	12	!! stp_RK3 : RK3 Shallow Water Eq. time-stepping
[12983]	13	!!----------------------------------------------------------------------
[13604]	14	USE stp_oce ! modules used in nemo_init and stp_RK3
[12983]	15	!
[13604]	16	USE domqco ! quasi-eulerian coordinate
	17	USE phycst ! physical constants
	18	USE usrdef_nam ! user defined namelist parameters
[12983]	19
	20	IMPLICIT NONE
	21	PRIVATE
	22
	23	PUBLIC stp_RK3 ! called by nemogcm.F90
[13604]	24
	25	! ! time level indices !
	26	INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs !: used by nemo_init
[12983]	27
	28	!! * Substitutions
	29	# include "do_loop_substitute.h90"
	30	# include "domzgr_substitute.h90"
	31	!!----------------------------------------------------------------------
	32	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
	33	!! $Id: step.F90 12614 2020-03-26 14:59:52Z gm $
	34	!! Software governed by the CeCILL license (see ./LICENSE)
	35	!!----------------------------------------------------------------------
	36	CONTAINS
	37
	38	SUBROUTINE stp_RK3( kstp )
	39	!!----------------------------------------------------------------------
	40	!! * ROUTINE stp_RK3 *
	41	!!
[13604]	42	!! ** Purpose : - RK3 Time stepping scheme for shallow water Eq.
[12983]	43	!!
[13604]	44	!! ** Method : 3rd order time stepping scheme which has 3 stages
	45	!! * Update calendar and forcings
	46	!! stage 1 : n ==> n+1/3 using variables at n
	47	!! - Compute the rhs of momentum
	48	!! - Time step ssh at Naa (n+1/3)
	49	!! - Time step u,v at Naa (n+1/3)
	50	!! - Swap time indices
	51	!! stage 2 : n ==> n+1/2 using variables at n and n+1/3
	52	!! - Compute the rhs of momentum
	53	!! - Time step ssh at Naa (n+1/2)
	54	!! - Time step u,v at Naa (n+1/2)
	55	!! - Swap time indices
	56	!! stage 3 : n ==> n+1 using variables at n and n+1/2
	57	!! - Compute the rhs of momentum
	58	!! - Time step ssh at Naa (n+1)
	59	!! - Time step u,v at Naa (n+1)
	60	!! - Swap time indices
	61	!! * Outputs and diagnostics
	62	!!
	63	!! NB: in stages 1 and 2 lateral mixing and forcing are not taken
	64	!! into account in the momentum RHS execpt if key_RK3all is used
[12983]	65	!!----------------------------------------------------------------------
[13604]	66	INTEGER, INTENT(in ) :: kstp ! ocean time-step index
	67	!
[12983]	68	INTEGER :: ji, jj, jk ! dummy loop indice
	69	INTEGER :: indic ! error indicator if < 0
	70	REAL(wp):: z1_2rho0, z5_6, z3_4 ! local scalars
[13769]	71	REAL(wp):: zue3a, zue3b, zua, zrhs_u ! local scalars
	72	REAL(wp):: zve3a, zve3b, zva, zrhs_v ! - -
[12983]	73	!! ---------------------------------------------------------------------
	74	!
	75	IF( ln_timing ) CALL timing_start('stp_RK3')
	76	!
	77	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	78	! model timestep
	79	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	80	!
	81	IF ( kstp == nit000 ) ww(:,:,:) = 0._wp ! initialize vertical velocity one for all to zero
	82
	83	!
	84	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	85	! update I/O and calendar
	86	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	87	indic = 0 ! reset to no error condition
	88
[13604]	89	IF( kstp == nit000 ) THEN ! initialize IOM context
	90	CALL iom_init( cxios_context, ld_closedef=.FALSE. ) ! for model grid (including possible AGRIF zoom)
[12983]	91	CALL iom_init_closedef
	92	ENDIF
	93	IF( kstp /= nit000 ) CALL day( kstp ) ! Calendar (day was already called at nit000 in day_init)
	94	CALL iom_setkt( kstp - nit000 + 1, cxios_context ) ! tell IOM we are at time step kstp
	95
	96	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
[13604]	97	! Update external forcing (SWE: surface boundary condition only)
[12983]	98	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	99
[13604]	100	CALL sbc ( kstp, Nbb, Nnn ) ! Sea Boundary Condition
	101
[12983]	102	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
[13604]	103	! Ocean physics update (SWE: eddy viscosity only)
[12983]	104	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[13604]	105
[12983]	106	IF( l_ldfdyn_time ) CALL ldf_dyn( kstp, Nbb ) ! eddy viscosity coeff.
	107
	108	!======================================================================
	109	!======================================================================
	110	! ===== RK3 =====
	111	!======================================================================
	112	!======================================================================
	113
	114
	115	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
[13604]	116	! RK3 1st stage Ocean dynamics : u, v, ssh
[12983]	117	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	118	rDt = rn_Dt / 3._wp
	119	r1_Dt = 1._wp / rDt
[13604]	120	!
	121	! !== RHS of the momentum Eq. ==!
	122	!
	123	uu(:,:,:,Nrhs) = 0._wp ! set dynamics trends to zero
	124	vv(:,:,:,Nrhs) = 0._wp
[12983]	125
[13604]	126	CALL dyn_adv( kstp, Nbb, Nbb, uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS
	127	CALL dyn_vor( kstp, Nbb, uu, vv, Nrhs ) ! vorticity ==> RHS
[12983]	128	#if defined key_RK3all
[13604]	129	CALL dyn_ldf( kstp, Nbb, Nbb, uu, vv, Nrhs ) ! lateral mixing
[12983]	130	#endif
[13769]	131	z5_6 = 5._wp/6._wp
	132	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[13604]	133	! ! horizontal pressure gradient
[13769]	134	zrhs_u = - grav * ( ssh(ji+1,jj,Nbb) - ssh(ji,jj,Nbb) ) * r1_e1u(ji,jj)
	135	zrhs_v = - grav * ( ssh(ji,jj+1,Nbb) - ssh(ji,jj,Nbb) ) * r1_e2v(ji,jj)
	136	#if defined key_RK3all
	137	! ! wind stress and layer friction
	138	zrhs_u = zrhs_u + r1_rho0 * ( z5_6utau_b(ji,jj) + (1._wp - z5_6)utau(ji,jj) ) / e3u(ji,jj,jk,Nbb) &
	139	& - rn_rfr * uu(ji,jj,jk,Nbb)
	140	zrhs_v = zrhs_v + r1_rho0 * ( z5_6vtau_b(ji,jj) + (1._wp - z5_6)vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb) &
	141	& - rn_rfr * vv(ji,jj,jk,Nbb)
	142	#endif
	143	! ! ==> RHS
	144	uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
	145	vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
[12983]	146	END_3D
	147	!
[13604]	148	! !== Time stepping of ssh Eq. ==! (and update r3_Naa)
	149	!
	150	CALL ssh_nxt( kstp, Nbb, Nbb, ssh, Naa ) ! after ssh
[13769]	151	! ! after ssh/h_0 ratio
[13604]	152	CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
	153	!
	154	! !== Time stepping of momentum Eq. ==!
	155	!
	156	IF( ln_dynadv_vec ) THEN ! vector invariant form : applied on velocity
[13769]	157	DO_3D( 0, 0, 0, 0, 1,jpkm1)
[12983]	158	uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
	159	vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
	160	END_3D
	161	ELSE
[13769]	162	DO_3D( 0, 0, 0, 0, 1,jpkm1) ! flux form : applied on thickness weighted velocity
[13604]	163	zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
	164	zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
	165	zue3a = zue3b + rDt * e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
	166	zve3a = zve3b + rDt * e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
	167	!
	168	uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)
	169	vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
[12983]	170	END_3D
	171	ENDIF
[13604]	172	!
[13769]	173	CALL lbc_lnk_multi( 'stp_RK3', uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. )
	174	!
[13604]	175	! !== Swap time levels ==!
[12983]	176	Nrhs= Nnn
	177	Nnn = Naa
	178	Naa = Nrhs
	179
	180	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	181	! RK3 2nd stage Ocean dynamics : hdiv, ssh, e3, u, v, w
	182	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	183	rDt = rn_Dt / 2._wp
	184	r1_Dt = 1._wp / rDt
[13604]	185	!
	186	! !== RHS of the momentum Eq. ==!
	187	!
	188	uu(:,:,:,Nrhs) = 0._wp ! set dynamics trends to zero
	189	vv(:,:,:,Nrhs) = 0._wp
	190	CALL dyn_adv( kstp, Nbb, Nnn, uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS
	191	CALL dyn_vor( kstp, Nnn, uu, vv, Nrhs ) ! vorticity ==> RHS
[12983]	192	#if defined key_RK3all
[13604]	193	CALL dyn_ldf( kstp, Nbb, Nbb, uu, vv, Nrhs ) ! lateral mixing
[12983]	194	#endif
	195	!
[13769]	196	z3_4 = 3._wp/4._wp
[13295]	197	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[13604]	198	! ! horizontal pressure gradient
[13769]	199	zrhs_u = - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
	200	zrhs_v = - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
[12983]	201	#if defined key_RK3all
[13769]	202	! ! wind stress and layer friction
	203	zrhs_u = zrhs_u + r1_rho0 * ( z3_4utau_b(ji,jj) + (1._wp - z3_4)utau(ji,jj) ) / e3u(ji,jj,jk,Nnn) &
	204	& - rn_rfr * uu(ji,jj,jk,Nbb)
	205	zrhs_v = zrhs_v + r1_rho0 * ( z3_4vtau_b(ji,jj) + (1._wp - z3_4)vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn) &
	206	& - rn_rfr * vv(ji,jj,jk,Nbb)
	207	#endif
	208	! ! ==> RHS
	209	uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
	210	vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
[12983]	211	END_3D
[13604]	212	!
	213	! !== Time stepping of ssh Eq. ==! (and update r3_Naa)
	214	!
	215	CALL ssh_nxt( kstp, Nbb, Nnn, ssh, Naa ) ! after ssh
[13769]	216	! ! after ssh/h_0 ratio
[13604]	217	CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
	218	!
	219	! !== Time stepping of momentum Eq. ==!
	220	!
	221	IF( ln_dynadv_vec ) THEN ! vector invariant form : applied on velocity
[13769]	222	DO_3D( 0, 0, 0, 0, 1,jpkm1)
[12983]	223	uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
	224	vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
	225	END_3D
	226	ELSE
[13769]	227	DO_3D( 0, 0, 0, 0, 1,jpkm1) ! flux form : applied on thickness weighted velocity
[13604]	228	zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
	229	zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
	230	zue3a = zue3b + rDt * e3u(ji,jj,jk,Nnn) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
	231	zve3a = zve3b + rDt * e3v(ji,jj,jk,Nnn) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
	232	!
	233	uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)
	234	vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
[12983]	235	END_3D
	236	ENDIF
[13604]	237	!
[13769]	238	CALL lbc_lnk_multi( 'stp_RK3', uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. )
	239	!
[13604]	240	! !== Swap time levels ==!
[12983]	241	Nrhs= Nnn
	242	Nnn = Naa
	243	Naa = Nrhs
	244
	245	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	246	! RK3 3rd stage Ocean dynamics : hdiv, ssh, e3, u, v, w
	247	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
	248	rDt = rn_Dt
	249	r1_Dt = 1._wp / rDt
[13604]	250	!
	251	! !== RHS of the momentum Eq. ==!
	252	!
	253	uu(:,:,:,Nrhs) = 0._wp ! set dynamics trends to zero
	254	vv(:,:,:,Nrhs) = 0._wp
	255	!
	256	CALL dyn_adv( kstp, Nbb, Nnn, uu, vv, Nrhs ) ! advection (VF or FF) ==> RHS
	257	CALL dyn_vor( kstp, Nnn, uu, vv, Nrhs ) ! vorticity ==> RHS
	258	CALL dyn_ldf( kstp, Nbb, Nnn, uu, vv, Nrhs ) ! lateral mixing
[12983]	259
[13769]	260	z1_2rho0 = 0.5_wp * r1_rho0
	261	DO_3D( 0, 0, 0, 0, 1,jpkm1 )
[13604]	262	! ! horizontal pressure gradient
[13769]	263	zrhs_u = - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
	264	zrhs_v = - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
	265	! ! wind stress and layer friction
	266	zrhs_u = zrhs_u + z1_2rho0 * ( utau_b(ji,jj) + utau(ji,jj) ) / e3u(ji,jj,jk,Nnn) &
	267	& - rn_rfr * uu(ji,jj,jk,Nbb)
	268	zrhs_v = zrhs_v + z1_2rho0 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / e3v(ji,jj,jk,Nnn) &
	269	& - rn_rfr * vv(ji,jj,jk,Nbb)
	270	! ! ==> RHS
	271	uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + zrhs_u
	272	vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
[12983]	273	END_3D
[13604]	274	!
	275	! !== Time stepping of ssh Eq. ==! (and update r3_Naa)
	276	!
	277	CALL ssh_nxt( kstp, Nbb, Nnn, ssh, Naa ) ! after ssh
[13769]	278	! ! after ssh/h_0 ratio
[13604]	279	CALL dom_qco_r3c( ssh(:,:,Naa), r3t(:,:,Naa), r3u(:,:,Naa), r3v(:,:,Naa), r3f(:,:) )
	280	!
	281	! !== Time stepping of momentum Eq. ==!
	282	!
	283	IF( ln_dynadv_vec ) THEN ! vector invariant form : applied on velocity
[13769]	284	DO_3D( 0, 0, 0, 0, 1,jpkm1)
	285	uu(ji,jj,jk,Naa) = uu(ji,jj,jk,Nbb) + rDt * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
	286	vv(ji,jj,jk,Naa) = vv(ji,jj,jk,Nbb) + rDt * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
[12983]	287	END_3D
	288	!
[13604]	289	ELSE ! flux form : applied on thickness weighted velocity
[13769]	290	DO_3D( 0, 0, 0, 0, 1,jpkm1)
[12983]	291	zue3b = e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nbb)
	292	zve3b = e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nbb)
[13604]	293	zue3a = zue3b + rDt * e3u(ji,jj,jk,Nbb) * uu(ji,jj,jk,Nrhs) * umask(ji,jj,jk)
	294	zve3a = zve3b + rDt * e3v(ji,jj,jk,Nbb) * vv(ji,jj,jk,Nrhs) * vmask(ji,jj,jk)
[12983]	295	!
[13604]	296	uu(ji,jj,jk,Naa) = zue3a / e3u(ji,jj,jk,Naa)
	297	vv(ji,jj,jk,Naa) = zve3a / e3v(ji,jj,jk,Naa)
[12983]	298	END_3D
	299	ENDIF
[13604]	300	!
[13769]	301	CALL lbc_lnk_multi( 'stp_RK3', uu(:,:,:,Naa), 'U', -1., vv(:,:,:,Naa), 'V', -1. )
	302	!
[13604]	303	! !== Swap time levels ==!
	304	!
[12983]	305	Nrhs = Nbb
	306	Nbb = Naa
	307	Naa = Nrhs
[13604]	308
[12983]	309	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
[14179]	310	! diagnostics and outputs at Nbb (i.e. the just computed time step)
[12983]	311	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[13604]	312
[14179]	313	IF( ln_diacfl ) CALL dia_cfl ( kstp, Nbb ) ! Courant number diagnostics
	314	CALL dia_wri ( kstp, Nbb ) ! ocean model: outputs
[12983]	315	!
[14179]	316	IF( lrst_oce ) CALL rst_write ( kstp, Nbb, Nbb ) ! write output ocean restart file
[12983]	317
	318	!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
	319	! Control
	320	!<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
[14179]	321	CALL stp_ctl_SWE ( kstp , Nbb )
[13604]	322
[12983]	323	IF( kstp == nit000 ) THEN ! 1st time step only
	324	CALL iom_close( numror ) ! close input ocean restart file
	325	IF(lwm) CALL FLUSH ( numond ) ! flush output namelist oce
	326	IF(lwm .AND. numoni /= -1 ) CALL FLUSH ( numoni ) ! flush output namelist ice (if exist)
	327	ENDIF
	328
	329	!
[14239]	330	#if defined key_xios
[12983]	331	IF( kstp == nitend .OR. indic < 0 ) THEN
[13604]	332	CALL iom_context_finalize( cxios_context )
[12983]	333	ENDIF
	334	#endif
	335	!
	336	IF( ln_timing ) CALL timing_stop('stp_RK3')
	337	!
	338	END SUBROUTINE stp_RK3
[13604]	339
[12983]	340	!!======================================================================
[13769]	341	END MODULE stprk3

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

Download in other formats: