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.

bdytides.F90 in trunk/NEMO/OPA_SRC/BDY – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/BDY/bdytides.F90 @ 1473

Last change on this file since 1473 was 1462, checked in by rblod, 15 years ago
Bugs in bdytides, see ticket #441
Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 29.2 KB

Line
1	MODULE bdytides
2	!!======================================================================
3	!! * MODULE bdytides *
4	!! Ocean dynamics: Tidal forcing at open boundaries
5	!!======================================================================
6	!! History : 2.0 ! 2007-01 (D.Storkey) Original code
7	!! 2.3 ! 2008-01 (J.Holt) Add date correction. Origins POLCOMS v6.3 2007
8	!! 3.0 ! 2008-04 (NEMO team) add in the reference version
9	!!----------------------------------------------------------------------
10	#if defined key_bdy
11	!!----------------------------------------------------------------------
12	!! 'key_bdy' Unstructured Open Boundary Condition
13	!!----------------------------------------------------------------------
14	!! PUBLIC
15	!! tide_init : read of namelist
16	!! tide_data : read in and initialisation of tidal constituents at boundary
17	!! tide_update : calculation of tidal forcing at each timestep
18	!! PRIVATE
19	!! uvset :\
20	!! vday : \| Routines to correct tidal harmonics forcing for
21	!! shpen : \| start time of integration
22	!! ufset : \|
23	!! vset :/
24	!!----------------------------------------------------------------------
25	USE oce ! ocean dynamics and tracers
26	USE dom_oce ! ocean space and time domain
27	USE iom
28	USE in_out_manager ! I/O units
29	USE phycst ! physical constants
30	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
31	USE bdy_par ! Unstructured boundary parameters
32	USE bdy_oce ! ocean open boundary conditions
33	USE daymod ! calendar
34
35	IMPLICIT NONE
36	PRIVATE
37
38	PUBLIC tide_init ! routine called in bdyini
39	PUBLIC tide_data ! routine called in bdyini
40	PUBLIC tide_update ! routine called in bdydyn
41
42	LOGICAL, PUBLIC :: ln_tide_date !: =T correct tide phases and amplitude for model start date
43
44	INTEGER, PARAMETER :: jptides_max = 15 !: Max number of tidal contituents
45	INTEGER :: ntide !: Actual number of tidal constituents
46
47	CHARACTER(len=80), PUBLIC :: filtide !: Filename root for tidal input files
48	CHARACTER(len= 4), PUBLIC, DIMENSION(jptides_max) :: tide_cpt !: Names of tidal components used.
49
50	INTEGER , DIMENSION(jptides_max) :: nindx !: ???
51	REAL(wp), DIMENSION(jptides_max) :: tide_speed !: Phase speed of tidal constituent (deg/hr)
52
53	REAL(wp), DIMENSION(jpbdim,jptides_max) :: ssh1, ssh2 !: Tidal constituents : SSH
54	REAL(wp), DIMENSION(jpbdim,jptides_max) :: u1 , u2 !: Tidal constituents : U
55	REAL(wp), DIMENSION(jpbdim,jptides_max) :: v1 , v2 !: Tidal constituents : V
56
57	!!----------------------------------------------------------------------
58	!! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
59	!! $Id$
60	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
61	!!----------------------------------------------------------------------
62
63	CONTAINS
64
65	SUBROUTINE tide_init
66	!!----------------------------------------------------------------------
67	!! * SUBROUTINE tide_init *
68	!!
69	!! ** Purpose : - Read in namelist for tides
70	!!
71	!!----------------------------------------------------------------------
72	INTEGER :: itide ! dummy loop index
73	!!
74	NAMELIST/namtide/ln_tide_date, filtide, tide_cpt, tide_speed
75	!!----------------------------------------------------------------------
76
77	IF(lwp) WRITE(numout,*)
78	IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries'
79	IF(lwp) WRITE(numout,*) '~~~~~~~~~'
80
81	! Namelist namtide : tidal harmonic forcing at open boundaries
82	ln_tide_date = .false.
83	filtide(:) = ''
84	tide_speed(:) = 0.0
85	tide_cpt(:) = ''
86	REWIND( numnam ) ! Read namelist parameters
87	READ ( numnam, namtide )
88	! ! Count number of components specified
89	ntide = jptides_max
90	itide = 1
91	DO WHILE( tide_cpt(itide) /= '' )
92	ntide = itide
93	itide = itide + 1
94	END DO
95	! ! find constituents in standard list
96	DO itide = 1, ntide
97	nindx(itide) = 0
98	IF( TRIM( tide_cpt(itide) ) == 'Q1' ) nindx(itide) = 1
99	IF( TRIM( tide_cpt(itide) ) == 'O1' ) nindx(itide) = 2
100	IF( TRIM( tide_cpt(itide) ) == 'P1' ) nindx(itide) = 3
101	IF( TRIM( tide_cpt(itide) ) == 'S1' ) nindx(itide) = 4
102	IF( TRIM( tide_cpt(itide) ) == 'K1' ) nindx(itide) = 5
103	IF( TRIM( tide_cpt(itide) ) == '2N2' ) nindx(itide) = 6
104	IF( TRIM( tide_cpt(itide) ) == 'MU2' ) nindx(itide) = 7
105	IF( TRIM( tide_cpt(itide) ) == 'N2' ) nindx(itide) = 8
106	IF( TRIM( tide_cpt(itide) ) == 'NU2' ) nindx(itide) = 9
107	IF( TRIM( tide_cpt(itide) ) == 'M2' ) nindx(itide) = 10
108	IF( TRIM( tide_cpt(itide) ) == 'L2' ) nindx(itide) = 11
109	IF( TRIM( tide_cpt(itide) ) == 'T2' ) nindx(itide) = 12
110	IF( TRIM( tide_cpt(itide) ) == 'S2' ) nindx(itide) = 13
111	IF( TRIM( tide_cpt(itide) ) == 'K2' ) nindx(itide) = 14
112	IF( TRIM( tide_cpt(itide) ) == 'M4' ) nindx(itide) = 15
113	IF( nindx(itide) == 0 .AND. lwp ) THEN
114	WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list'
115	CALL ctl_warn( ctmp1 )
116	ENDIF
117	END DO
118	! ! Parameter control and print
119	IF( ntide < 1 ) THEN
120	CALL ctl_stop( ' Did not find any tidal components in namelist namtide' )
121	ELSE
122	IF(lwp) WRITE(numout,*) ' Namelist namtide : tidal harmonic forcing at open boundaries'
123	IF(lwp) WRITE(numout,*) ' tidal components specified ', ntide
124	IF(lwp) WRITE(numout,*) ' ', tide_cpt(1:ntide)
125	IF(lwp) WRITE(numout,*) ' associated phase speeds (deg/hr) : '
126	IF(lwp) WRITE(numout,*) ' ', tide_speed(1:ntide)
127	ENDIF
128
129	! Initialisation of tidal harmonics arrays
130	sshtide(:) = 0.e0
131	utide (:) = 0.e0
132	vtide (:) = 0.e0
133	!
134	END SUBROUTINE tide_init
135
136
137	SUBROUTINE tide_data
138	!!----------------------------------------------------------------------
139	!! * SUBROUTINE tide_data *
140	!!
141	!! ** Purpose : - Read in tidal harmonics data and adjust for the start
142	!! time of the model run.
143	!!
144	!!----------------------------------------------------------------------
145	INTEGER :: itide, igrd, ib ! dummy loop indices
146	CHARACTER(len=80) :: clfile ! full file name for tidal input file
147	INTEGER :: ipi, ipj, inum, idvar ! temporary integers (netcdf read)
148	INTEGER, DIMENSION(3) :: lendta=0 ! length of data in the file (note may be different from nblendta!)
149	REAL(wp) :: z_arg, z_atde, z_btde, z1t, z2t
150	REAL(wp), DIMENSION(jpbdta,1) :: zdta ! temporary array for data fields
151	REAL(wp), DIMENSION(jptides_max) :: z_vplu, z_ftc
152	!!------------------------------------------------------------------------------
153
154	! Open files and read in tidal forcing data
155	! -----------------------------------------
156
157	ipj = 1
158
159	DO itide = 1, ntide
160	! ! SSH fields
161	clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc'
162	IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
163	CALL iom_open( clfile, inum )
164	igrd = 1
165	IF( nblendta(igrd) <= 0 ) THEN
166	idvar = iom_varid( inum,'z1' )
167	IF(lwp) WRITE(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar)
168	nblendta(igrd) = iom_file(inum)%dimsz(1,idvar)
169	WRITE(numout,*) 'Dim size for z1 is ', nblendta(igrd)
170	ENDIF
171	ipi = nblendta(igrd)
172	CALL iom_get( inum, jpdom_unknown, 'z1', zdta(1:ipi,1:ipj) )
173	DO ib = 1, nblenrim(igrd)
174	ssh1(ib,itide) = zdta(nbmap(ib,igrd),1)
175	END DO
176	CALL iom_get( inum, jpdom_unknown, 'z2', zdta(1:ipi,1:ipj) )
177	DO ib = 1, nblenrim(igrd)
178	ssh2(ib,itide) = zdta(nbmap(ib,igrd),1)
179	END DO
180	CALL iom_close( inum )
181	!
182	! ! U fields
183	clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc'
184	IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
185	CALL iom_open( clfile, inum )
186	igrd = 2
187	IF( lendta(igrd) <= 0 ) THEN
188	idvar = iom_varid( inum,'u1' )
189	lendta(igrd) = iom_file(inum)%dimsz(1,idvar)
190	WRITE(numout,*) 'Dim size for u1 is ',lendta(igrd)
191	ENDIF
192	ipi = lendta(igrd)
193	CALL iom_get( inum, jpdom_unknown, 'u1', zdta(1:ipi,1:ipj) )
194	DO ib = 1, nblenrim(igrd)
195	u1(ib,itide) = zdta(nbmap(ib,igrd),1)
196	END DO
197	CALL iom_get( inum, jpdom_unknown, 'u2', zdta(1:ipi,1:ipj) )
198	DO ib = 1, nblenrim(igrd)
199	u2(ib,itide) = zdta(nbmap(ib,igrd),1)
200	END DO
201	CALL iom_close( inum )
202	!
203	! ! V fields
204	clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc'
205	if(lwp) write(numout,*) 'Reading data from file ', clfile
206	CALL iom_open( clfile, inum )
207	igrd = 3
208	IF( lendta(igrd) <= 0 ) THEN
209	idvar = iom_varid( inum,'v1' )
210	lendta(igrd) = iom_file(inum)%dimsz(1,idvar)
211	WRITE(numout,*) 'Dim size for v1 is ', lendta(igrd)
212	ENDIF
213	ipi = lendta(igrd)
214	CALL iom_get( inum, jpdom_unknown, 'v1', zdta(1:ipi,1:ipj) )
215	DO ib = 1, nblenrim(igrd)
216	v1(ib,itide) = zdta(nbmap(ib,igrd),1)
217	END DO
218	CALL iom_get( inum, jpdom_unknown, 'v2', zdta(1:ipi,1:ipj) )
219	DO ib=1, nblenrim(igrd)
220	v2(ib,itide) = zdta(nbmap(ib,igrd),1)
221	END DO
222	CALL iom_close( inum )
223	!
224	END DO ! end loop on tidal components
225
226	IF( ln_tide_date ) THEN ! correct for date factors
227
228	!! used nmonth, nyear and nday from daymod....
229	! Calculate date corrects for 15 standard consituents
230	! This is the initialisation step, so nday, nmonth, nyear are the
231	! initial date/time of the integration.
232	print *, nday,nmonth,nyear
233	nyear = int(ndate0 / 10000 ) ! initial year
234	nmonth = int((ndate0 - nyear * 10000 ) / 100 ) ! initial month
235	nday = int(ndate0 - nyear * 10000 - nmonth * 100)
236
237	CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu )
238
239	IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear
240
241	DO itide = 1, ntide ! loop on tidal components
242	!
243	IF( nindx(itide) /= 0 ) THEN
244	!!gm use rpi and rad global variable
245	z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0
246	z_atde=z_ftc(nindx(itide))*cos(z_arg)
247	z_btde=z_ftc(nindx(itide))*sin(z_arg)
248	IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), tide_speed(itide), &
249	& z_ftc(nindx(itide)), z_vplu(nindx(itide))
250	ELSE
251	z_atde = 1.0_wp
252	z_btde = 0.0_wp
253	ENDIF
254	! ! elevation
255	igrd = 1
256	DO ib = 1, nblenrim(igrd)
257	z1t = z_atde * ssh1(ib,itide) + z_btde * ssh2(ib,itide)
258	z2t = z_atde * ssh2(ib,itide) - z_btde * ssh1(ib,itide)
259	ssh1(ib,itide) = z1t
260	ssh2(ib,itide) = z2t
261	END DO
262	! ! u
263	igrd = 2
264	DO ib = 1, nblenrim(igrd)
265	z1t = z_atde * u1(ib,itide) + z_btde * u2(ib,itide)
266	z2t = z_atde * u2(ib,itide) - z_btde * u1(ib,itide)
267	u1(ib,itide) = z1t
268	u2(ib,itide) = z2t
269	END DO
270	! ! v
271	igrd = 3
272	DO ib = 1, nblenrim(igrd)
273	z1t = z_atde * v1(ib,itide) + z_btde * v2(ib,itide)
274	z2t = z_atde * v2(ib,itide) - z_btde * v1(ib,itide)
275	v1(ib,itide) = z1t
276	v2(ib,itide) = z2t
277	END DO
278	!
279	END DO ! end loop on tidal components
280	!
281	ENDIF ! date correction
282	!
283	END SUBROUTINE tide_data
284
285
286	SUBROUTINE tide_update ( kt, jit )
287	!!----------------------------------------------------------------------
288	!! * SUBROUTINE tide_update *
289	!!
290	!! ** Purpose : - Add tidal forcing to sshbdy, ubtbdy and vbtbdy arrays.
291	!!
292	!!----------------------------------------------------------------------
293	INTEGER, INTENT( in ) :: kt ! Main timestep counter
294	!!gm doctor jit ==> kit
295	INTEGER, INTENT( in ) :: jit ! Barotropic timestep counter (for timesplitting option)
296	!!
297	INTEGER :: itide, igrd, ib ! dummy loop indices
298	REAL(wp) :: z_arg, z_sarg !
299	REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost
300	!!----------------------------------------------------------------------
301
302	! Note tide phase speeds are in deg/hour, so we need to convert the
303	! elapsed time in seconds to hours by dividing by 3600.0
304	IF( jit == 0 ) THEN
305	z_arg = kt * rdt * rad / 3600.0
306	ELSE ! we are in a barotropic subcycle (for timesplitting option)
307	! z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,lwp) ) * rad / 3600.0
308	z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,wp) ) * rad / 3600.0
309	ENDIF
310
311	DO itide = 1, ntide
312	z_sarg = z_arg * tide_speed(itide)
313	z_cost(itide) = COS( z_sarg )
314	z_sist(itide) = SIN( z_sarg )
315	END DO
316
317	! summing of tidal constituents into BDY arrays
318	sshtide(:) = 0.0
319	utide (:) = 0.0
320	vtide (:) = 0.0
321	!
322	DO itide = 1, ntide
323	igrd=1 ! SSH on tracer grid.
324	DO ib = 1, nblenrim(igrd)
325	sshtide(ib) =sshtide(ib)+ ssh1(ib,itide)z_cost(itide) + ssh2(ib,itide)z_sist(itide)
326	! if(lwp) write(numout,*) 'z',ib,itide,sshtide(ib), ssh1(ib,itide),ssh2(ib,itide)
327	END DO
328	igrd=2 ! U grid
329	DO ib=1, nblenrim(igrd)
330	utide(ib) = utide(ib)+ u1(ib,itide)z_cost(itide) + u2(ib,itide)z_sist(itide)
331	! if(lwp) write(numout,*) 'u',ib,itide,utide(ib), u1(ib,itide),u2(ib,itide)
332	END DO
333	igrd=3 ! V grid
334	DO ib=1, nblenrim(igrd)
335	vtide(ib) = vtide(ib)+ v1(ib,itide)z_cost(itide) + v2(ib,itide)z_sist(itide)
336	! if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), v1(ib,itide),v2(ib,itide)
337	END DO
338	END DO
339	!
340	END SUBROUTINE tide_update
341
342	!!gm doctor naming of dummy argument variables!!! and all local variables
343
344	SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu )
345	!!----------------------------------------------------------------------
346	!! * SUBROUTINE uvset *
347	!!
348	!! ** Purpose : - adjust tidal forcing for date factors
349	!!
350	!!----------------------------------------------------------------------
351	implicit none
352	INTEGER, INTENT( in ) :: ihs ! Start time hours
353	INTEGER, INTENT( in ) :: iday ! start time days
354	INTEGER, INTENT( in ) :: imnth ! start time month
355	INTEGER, INTENT( in ) :: iyr ! start time year
356	!!
357	!!gm nc is jptides_max ????
358	INTEGER , PARAMETER :: nc =15 ! maximum number of constituents
359	CHARACTER(len=8), DIMENSION(nc) :: cname
360	INTEGER :: year, vd, ivdy, ndc, i, k
361	REAL(wp) :: ss, h, p, en, p1, rtd
362	REAL(wp), DIMENSION(nc) :: f ! nodal correction
363	REAL(wp), DIMENSION(nc) :: z_vplu ! phase correction
364	REAL(wp), DIMENSION(nc) :: u, v, zig
365	!!
366	DATA cname/ 'q1' , 'o1' , 'p1' , 's1' , 'k1' , &
367	& '2n2' , 'mu2' , 'n2' , 'nu2' , 'm2' , &
368	& 'l2' , 't2' , 's2' , 'k2' , 'm4' /
369	DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878, .2625161701, &
370	& .4868657873, .4881373225, .4963669182, .4976384533, .5058680490, &
371	& .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098 /
372	!!----------------------------------------------------------------------
373	!
374	! ihs - start time gmt on ...
375	! iday/imnth/iyr - date e.g. 12/10/87
376	!
377	CALL vday(iday,imnth,iyr,ivdy)
378	vd = ivdy
379	!
380	!rp note change of year number for d. blackman shpen
381	!rp if(iyr.ge.1000) year=iyr-1900
382	!rp if(iyr.lt.1000) year=iyr
383	year = iyr
384	!
385	!.....year = year of required data
386	!.....vd = day of required data..set up for 0000gmt day year
387	ndc = nc
388	!.....ndc = number of constituents allowed
389	!!gm use rpi ?
390	rtd = 360.0 / 6.2831852
391	DO i = 1, ndc
392	zig(i) = zig(i)*rtd
393	! sigo(i)= zig(i)
394	END DO
395
396	!!gm try to avoid RETURN in F90
397	IF( year == 0 ) RETURN
398
399	CALL shpen( year, vd, ss, h , p , en, p1 )
400	CALL ufset( p , en, u , f )
401	CALL vset ( ss , h , p , en, p1, v )
402	!
403	DO k = 1, nc
404	z_vplu(k) = v(k) + u(k)
405	z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k)
406	DO WHILE( z_vplu(k) < 0 )
407	z_vplu(k) = z_vplu(k) + 360.0
408	END DO
409	DO WHILE( z_vplu(k) > 360. )
410	z_vplu(k) = z_vplu(k) - 360.0
411	END DO
412	END DO
413	!
414	END SUBROUTINE uvset
415
416
417	SUBROUTINE vday( iday, imnth, iy, ivdy )
418	!!----------------------------------------------------------------------
419	!! * SUBROUTINE vday *
420	!!
421	!! ** Purpose : - adjust tidal forcing for date factors
422	!!
423	!!----------------------------------------------------------------------
424	INTEGER, INTENT(in ) :: iday, imnth, iy ! ????
425	INTEGER, INTENT( out) :: ivdy ! ???
426	!!
427	INTEGER :: iyr
428	!!------------------------------------------------------------------------------
429
430	!!gm nday_year in day mode is the variable compiuted here, no?
431	!!gm nday_year , & !: curent day counted from jan 1st of the current year
432
433	!calculate day number in year from day/month/year
434	if(imnth.eq.1) ivdy=iday
435	if(imnth.eq.2) ivdy=iday+31
436	if(imnth.eq.3) ivdy=iday+59
437	if(imnth.eq.4) ivdy=iday+90
438	if(imnth.eq.5) ivdy=iday+120
439	if(imnth.eq.6) ivdy=iday+151
440	if(imnth.eq.7) ivdy=iday+181
441	if(imnth.eq.8) ivdy=iday+212
442	if(imnth.eq.9) ivdy=iday+243
443	if(imnth.eq.10) ivdy=iday+273
444	if(imnth.eq.11) ivdy=iday+304
445	if(imnth.eq.12) ivdy=iday+334
446	iyr=iy
447	if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1
448	if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1
449	if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1
450	!
451	END SUBROUTINE vday
452
453	!!doctor norme for dummy arguments
454
455	SUBROUTINE shpen( year, vd, s, h, p, en, p1 )
456	!!----------------------------------------------------------------------
457	!! * SUBROUTINE shpen *
458	!!
459	!! ** Purpose : - calculate astronomical arguments for tides
460	!! this version from d. blackman 30 nove 1990
461	!!
462	!!----------------------------------------------------------------------
463	!!gm add INTENT in, out or inout.... DOCTOR name....
464	!!gm please do not use variable name with a single letter: impossible to search in a code
465	INTEGER :: year,vd
466	REAL(wp) :: s,h,p,en,p1
467	!!
468	INTEGER :: yr,ilc,icent,it,iday,ild,ipos,nn,iyd
469	REAL(wp) :: cycle,t,td,delt(84),delta,deltat
470	!!
471	DATA delt /-5.04, -3.90, -2.87, -0.58, 0.71, 1.80, &
472	& 3.08, 4.63, 5.86, 7.21, 8.58, 10.50, 12.10, &
473	& 12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39, &
474	& 19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68, &
475	& 22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63, &
476	& 23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80, &
477	& 24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89, &
478	& 27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96, &
479	& 31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39, &
480	& 33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87, &
481	& 38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00, &
482	& 45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81, &
483	& 52.57 /
484	!!----------------------------------------------------------------------
485
486	cycle = 360.0
487	ilc = 0
488	icent = year / 100
489	yr = year - icent * 100
490	t = icent - 20
491	!
492	! for the following equations
493	! time origin is fixed at 00 hr of jan 1st,2000.
494	! see notes by cartwright
495	!
496	!!gm old coding style, use CASE instead and avoid GOTO (obsolescence in fortran 90)
497	!!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90
498	it = icent - 20
499	if (it) 1,2,2
500	1 iday = it/4 -it
501	go to 3
502	2 iday = (it+3)/4 - it
503	!
504	! t is in julian century
505	! correction in gegorian calander where only century year divisible
506	! by 4 is leap year.
507	!
508	3 continue
509	!
510	td = 0.0
511	!
512	!!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90
513	if (yr) 4,5,4
514	!
515	4 iyd = 365*yr
516	ild = (yr-1)/4
517	if((icent - (icent/4)*4) .eq. 0) ilc = 1
518	td = iyd + ild + ilc
519	!
520	5 td = td + iday + vd -1.0 - 0.5
521	t = t + (td/36525.0)
522	!
523	ipos=year-1899
524	if (ipos .lt. 0) go to 7
525	if (ipos .gt. 83) go to 6
526	!
527	delta = (delt(ipos+1)+delt(ipos))/2.0
528	go to 7
529	!
530	6 delta= (65.0-50.5)/20.0*(year-1980)+50.5
531	!
532	7 deltat = delta * 1.0e-6
533	!
534	!!gm precision of the computation : example for s it should be replace by:
535	!!gm s = 218.3165 + (481267.8813 - 0.0016t)t + 152.0*deltat ==> more precise modify the last digits results
536	s = 218.3165 + 481267.8813t - 0.0016tt + 152.0deltat
537	h = 280.4661 + 36000.7698 t + 0.0003tt + 11.0deltat
538	p = 83.3535 + 4069.0139 t - 0.0103t*t + deltat
539	en = 234.9555 + 1934.1363 t - 0.0021t*t + deltat
540	p1 = 282.9384 + 1.7195 t + 0.0005t*t
541	!
542	nn = s / cycle
543	s = s - nn * cycle
544	IF( s < 0.e0 ) s = s + cycle
545	!
546	nn = h / cycle
547	h = h - cycle * nn
548	IF( h < 0.e0 ) h = h + cycle
549	!
550	nn = p / cycle
551	p = p - cycle * nn
552	IF( p < 0.e0) p = p + cycle
553	!
554	nn = en / cycle
555	en = en - cycle * nn
556	IF( en < 0.e0 ) en = en + cycle
557	en = cycle - en
558	!
559	nn = p1 / cycle
560	p1 = p1 - nn * cycle
561	!
562	END SUBROUTINE shpen
563
564
565	SUBROUTINE ufset( p, cn, b, a )
566	!!----------------------------------------------------------------------
567	!! * SUBROUTINE ufset *
568	!!
569	!! ** Purpose : - calculate nodal parameters for the tides
570	!!
571	!!----------------------------------------------------------------------
572	!!gm doctor naming of dummy argument variables!!! and all local variables
573	!!gm nc is jptides_max ????
574	integer nc
575	parameter (nc=15)
576	REAL(wp) p,cn
577	!!
578
579	!!gm rad is already a public variable defined in phycst.F90 .... ==> doctor norme local real start with "z"
580	REAL(wp) :: w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad
581	REAL(wp) :: a(nc), b(nc)
582	INTEGER :: ndc, k
583	!!----------------------------------------------------------------------
584
585	ndc = nc
586
587	! a=f , b =u
588	! t is zero as compared to tifa.
589	!! use rad defined in phycst (i.e. add a USE phycst at the begining of the module
590	rad = 6.2831852d0/360.0
591	pw = p * rad
592	nw = cn * rad
593	w1 = cos( nw )
594	w2 = cos( 2*nw )
595	w3 = cos( 3*nw )
596	w4 = sin( nw )
597	w5 = sin( 2*nw )
598	w6 = sin( 3*nw )
599	w7 = 1. - 0.2505 * COS( 2pw ) - 0.1102 COS(2*pw-nw ) &
600	& - 0.156 * COS( 2pw-2nw ) - 0.037 * COS( nw )
601	w8 = - 0.2505 * SIN( 2pw ) - 0.1102 SIN(2*pw-nw ) &
602	& - 0.0156 * SIN( 2pw-2nw ) - 0.037 * SIN( nw )
603	!
604	a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3
605	b(1) = 0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6
606	! q1
607	a(2) = a(1)
608	b(2) = b(1)
609	! o1
610	a(3) = 1.0
611	b(3) = 0.0
612	! p1
613	a(4) = 1.0
614	b(4) = 0.0
615	! s1
616	a(5) = 1.0060+0.1150w1- 0.0088w2 +0.0006*w3
617	b(5) = -0.1546w4 + 0.0119w5 -0.0012*w6
618	! k1
619	a(6) =1.0004 -0.0373w1+ 0.0002w2
620	b(6) = -0.0374*w4
621	! 2n2
622	a(7) = a(6)
623	b(7) = b(6)
624	! mu2
625	a(8) = a(6)
626	b(8) = b(6)
627	! n2
628	a(9) = a(6)
629	b(9) = b(6)
630	! nu2
631	a(10) = a(6)
632	b(10) = b(6)
633	! m2
634	a(11) = SQRT( w7 * w7 + w8 * w8 )
635	b(11) = ATAN( w8 / w7 )
636	!!gmuse rpi instead of 3.141992 ??? true pi is rpi=3.141592653589793_wp ..... ????
637	IF( w7 < 0.e0 ) b(11) = b(11) + 3.141992
638	! l2
639	a(12) = 1.0
640	b(12) = 0.0
641	! t2
642	a(13)= a(12)
643	b(13)= b(12)
644	! s2
645	a(14) = 1.0241+0.2863w1+0.0083w2 -0.0015*w3
646	b(14) = -0.3096w4 + 0.0119w5 - 0.0007*w6
647	! k2
648	a(15) = a(6)*a(6)
649	b(15) = 2*b(6)
650	! m4
651	!!gm old coding, remove GOTO and label of lines
652	!!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90
653	DO 40 k = 1,ndc
654	b(k) = b(k)/rad
655	32 if (b(k)) 34,35,35
656	34 b(k) = b(k) + 360.0
657	go to 32
658	35 if (b(k)-360.0) 40,37,37
659	37 b(k) = b(k)-360.0
660	go to 35
661	40 continue
662	!
663	END SUBROUTINE ufset
664
665
666	SUBROUTINE vset( s,h,p,en,p1,v)
667	!!----------------------------------------------------------------------
668	!! * SUBROUTINE vset *
669	!!
670	!! ** Purpose : - calculate tidal phases for 0000gmt on start day of run
671	!!
672	!!----------------------------------------------------------------------
673	!!gm doctor naming of dummy argument variables!!! and all local variables
674	!!gm nc is jptides_max ????
675	!!gm en argument is not used: suppress it ?
676	integer nc
677	parameter (nc=15)
678	real(wp) s,h,p,en,p1
679	real(wp) v(nc)
680	!!
681	integer ndc, k
682	!!----------------------------------------------------------------------
683
684	ndc = nc
685	! v s are computed here.
686	v(1) =-3*s +h +p +270 ! Q1
687	v(2) =-2*s +h +270.0 ! O1
688	v(3) =-h +270 ! P1
689	v(4) =180 ! S1
690	v(5) =h +90.0 ! K1
691	v(6) =-4s +2h +2*p ! 2N2
692	v(7) =-4*(s-h) ! MU2
693	v(8) =-3s +2h +p ! N2
694	v(9) =-3s +4h -p ! MU2
695	v(10) =-2s +2h ! M2
696	v(11) =-s +2*h -p +180 ! L2
697	v(12) =-h +p1 ! T2
698	v(13) =0 ! S2
699	v(14) =h+h ! K2
700	v(15) =2*v(10) ! M4
701	!
702	!!gm old coding, remove GOTO and label of lines
703	!!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90
704	do 72 k = 1, ndc
705	69 if( v(k) ) 70,71,71
706	70 v(k) = v(k)+360.0
707	go to 69
708	71 if( v(k) - 360.0 ) 72,73,73
709	73 v(k) = v(k)-360.0
710	go to 71
711	72 continue
712	!
713	END SUBROUTINE vset
714
715	#else
716	!!----------------------------------------------------------------------
717	!! Dummy module NO Unstruct Open Boundary Conditions for tides
718	!!----------------------------------------------------------------------
719	!!gm are you sure we need to define filtide and tide_cpt ?
720	CHARACTER(len=80), PUBLIC :: filtide !: Filename root for tidal input files
721	CHARACTER(len=4 ), PUBLIC, DIMENSION(1) :: tide_cpt !: Names of tidal components used.
722
723	CONTAINS
724	SUBROUTINE tide_init ! Empty routine
725	END SUBROUTINE tide_init
726	SUBROUTINE tide_data ! Empty routine
727	END SUBROUTINE tide_data
728	SUBROUTINE tide_update( kt, kit ) ! Empty routine
729	WRITE(,) 'tide_update: You should not have seen this print! error?', kt, kit
730	END SUBROUTINE tide_update
731	#endif
732
733	!!======================================================================
734	END MODULE bdytides

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

Download in other formats: