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 branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Context Navigation

source: branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 3186

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

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

Download in other formats: