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_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 4474

Last change on this file since 4474 was 4436, checked in by trackstand2, 10 years ago
Added timing calls
Property svn:keywords set to `Id`
File size: 29.6 KB

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

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

Download in other formats: