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 @ 4409

Last change on this file since 4409 was 3211, checked in by spickles2, 13 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 29.5 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 )
[1125]	292	!!----------------------------------------------------------------------
	293	!! * SUBROUTINE tide_update *
	294	!!
[911]	295	!! ** Purpose : - Add tidal forcing to sshbdy, ubtbdy and vbtbdy arrays.
	296	!!
[1125]	297	!!----------------------------------------------------------------------
[911]	298	INTEGER, INTENT( in ) :: kt ! Main timestep counter
[1125]	299	!!gm doctor jit ==> kit
[911]	300	INTEGER, INTENT( in ) :: jit ! Barotropic timestep counter (for timesplitting option)
[1125]	301	!!
	302	INTEGER :: itide, igrd, ib ! dummy loop indices
	303	REAL(wp) :: z_arg, z_sarg !
	304	REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost
	305	!!----------------------------------------------------------------------
[911]	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
[1125]	309	IF( jit == 0 ) THEN
	310	z_arg = kt * rdt * rad / 3600.0
	311	ELSE ! we are in a barotropic subcycle (for timesplitting option)
[1462]	312	! z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,lwp) ) * rad / 3600.0
	313	z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,wp) ) * rad / 3600.0
[1125]	314	ENDIF
[911]	315
[1125]	316	DO itide = 1, ntide
	317	z_sarg = z_arg * tide_speed(itide)
	318	z_cost(itide) = COS( z_sarg )
	319	z_sist(itide) = SIN( z_sarg )
	320	END DO
[911]	321
[1125]	322	! summing of tidal constituents into BDY arrays
[911]	323	sshtide(:) = 0.0
[1125]	324	utide (:) = 0.0
	325	vtide (:) = 0.0
	326	!
	327	DO itide = 1, ntide
[2528]	328	igrd=4 ! SSH on tracer grid.
[1125]	329	DO ib = 1, nblenrim(igrd)
	330	sshtide(ib) =sshtide(ib)+ ssh1(ib,itide)z_cost(itide) + ssh2(ib,itide)z_sist(itide)
	331	! if(lwp) write(numout,*) 'z',ib,itide,sshtide(ib), ssh1(ib,itide),ssh2(ib,itide)
	332	END DO
[2528]	333	igrd=5 ! U grid
[1125]	334	DO ib=1, nblenrim(igrd)
	335	utide(ib) = utide(ib)+ u1(ib,itide)z_cost(itide) + u2(ib,itide)z_sist(itide)
	336	! if(lwp) write(numout,*) 'u',ib,itide,utide(ib), u1(ib,itide),u2(ib,itide)
	337	END DO
[2528]	338	igrd=6 ! V grid
[1125]	339	DO ib=1, nblenrim(igrd)
	340	vtide(ib) = vtide(ib)+ v1(ib,itide)z_cost(itide) + v2(ib,itide)z_sist(itide)
	341	! if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), v1(ib,itide),v2(ib,itide)
	342	END DO
	343	END DO
	344	!
	345	END SUBROUTINE tide_update
[911]	346
[1125]	347	!!gm doctor naming of dummy argument variables!!! and all local variables
[911]	348
[1125]	349	SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu )
	350	!!----------------------------------------------------------------------
	351	!! * SUBROUTINE uvset *
	352	!!
[911]	353	!! ** Purpose : - adjust tidal forcing for date factors
	354	!!
[1125]	355	!!----------------------------------------------------------------------
[911]	356	implicit none
	357	INTEGER, INTENT( in ) :: ihs ! Start time hours
	358	INTEGER, INTENT( in ) :: iday ! start time days
[1125]	359	INTEGER, INTENT( in ) :: imnth ! start time month
	360	INTEGER, INTENT( in ) :: iyr ! start time year
	361	!!
	362	!!gm nc is jptides_max ????
	363	INTEGER , PARAMETER :: nc =15 ! maximum number of constituents
[911]	364	CHARACTER(len=8), DIMENSION(nc) :: cname
[1125]	365	INTEGER :: year, vd, ivdy, ndc, i, k
	366	REAL(wp) :: ss, h, p, en, p1, rtd
	367	REAL(wp), DIMENSION(nc) :: f ! nodal correction
	368	REAL(wp), DIMENSION(nc) :: z_vplu ! phase correction
	369	REAL(wp), DIMENSION(nc) :: u, v, zig
	370	!!
	371	DATA cname/ 'q1' , 'o1' , 'p1' , 's1' , 'k1' , &
	372	& '2n2' , 'mu2' , 'n2' , 'nu2' , 'm2' , &
	373	& 'l2' , 't2' , 's2' , 'k2' , 'm4' /
	374	DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878, .2625161701, &
	375	& .4868657873, .4881373225, .4963669182, .4976384533, .5058680490, &
	376	& .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098 /
	377	!!----------------------------------------------------------------------
[911]	378	!
	379	! ihs - start time gmt on ...
	380	! iday/imnth/iyr - date e.g. 12/10/87
	381	!
[1125]	382	CALL vday(iday,imnth,iyr,ivdy)
	383	vd = ivdy
[911]	384	!
	385	!rp note change of year number for d. blackman shpen
	386	!rp if(iyr.ge.1000) year=iyr-1900
	387	!rp if(iyr.lt.1000) year=iyr
	388	year = iyr
	389	!
	390	!.....year = year of required data
	391	!.....vd = day of required data..set up for 0000gmt day year
	392	ndc = nc
	393	!.....ndc = number of constituents allowed
[1125]	394	!!gm use rpi ?
	395	rtd = 360.0 / 6.2831852
	396	DO i = 1, ndc
	397	zig(i) = zig(i)*rtd
	398	! sigo(i)= zig(i)
	399	END DO
	400
	401	!!gm try to avoid RETURN in F90
	402	IF( year == 0 ) RETURN
[911]	403
[1125]	404	CALL shpen( year, vd, ss, h , p , en, p1 )
	405	CALL ufset( p , en, u , f )
	406	CALL vset ( ss , h , p , en, p1, v )
	407	!
	408	DO k = 1, nc
	409	z_vplu(k) = v(k) + u(k)
	410	z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k)
	411	DO WHILE( z_vplu(k) < 0 )
	412	z_vplu(k) = z_vplu(k) + 360.0
	413	END DO
	414	DO WHILE( z_vplu(k) > 360. )
	415	z_vplu(k) = z_vplu(k) - 360.0
	416	END DO
	417	END DO
	418	!
[911]	419	END SUBROUTINE uvset
	420
	421
[1125]	422	SUBROUTINE vday( iday, imnth, iy, ivdy )
	423	!!----------------------------------------------------------------------
	424	!! * SUBROUTINE vday *
	425	!!
[911]	426	!! ** Purpose : - adjust tidal forcing for date factors
	427	!!
[1125]	428	!!----------------------------------------------------------------------
	429	INTEGER, INTENT(in ) :: iday, imnth, iy ! ????
	430	INTEGER, INTENT( out) :: ivdy ! ???
	431	!!
	432	INTEGER :: iyr
[911]	433	!!------------------------------------------------------------------------------
	434
[1125]	435	!!gm nday_year in day mode is the variable compiuted here, no?
	436	!!gm nday_year , & !: curent day counted from jan 1st of the current year
	437
	438	!calculate day number in year from day/month/year
[911]	439	if(imnth.eq.1) ivdy=iday
	440	if(imnth.eq.2) ivdy=iday+31
	441	if(imnth.eq.3) ivdy=iday+59
	442	if(imnth.eq.4) ivdy=iday+90
	443	if(imnth.eq.5) ivdy=iday+120
	444	if(imnth.eq.6) ivdy=iday+151
	445	if(imnth.eq.7) ivdy=iday+181
	446	if(imnth.eq.8) ivdy=iday+212
	447	if(imnth.eq.9) ivdy=iday+243
	448	if(imnth.eq.10) ivdy=iday+273
	449	if(imnth.eq.11) ivdy=iday+304
	450	if(imnth.eq.12) ivdy=iday+334
[1125]	451	iyr=iy
[911]	452	if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1
	453	if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1
	454	if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1
[1125]	455	!
	456	END SUBROUTINE vday
[911]	457
[1125]	458	!!doctor norme for dummy arguments
[911]	459
[1125]	460	SUBROUTINE shpen( year, vd, s, h, p, en, p1 )
	461	!!----------------------------------------------------------------------
	462	!! * SUBROUTINE shpen *
	463	!!
[911]	464	!! ** Purpose : - calculate astronomical arguments for tides
	465	!! this version from d. blackman 30 nove 1990
	466	!!
[1125]	467	!!----------------------------------------------------------------------
	468	!!gm add INTENT in, out or inout.... DOCTOR name....
	469	!!gm please do not use variable name with a single letter: impossible to search in a code
	470	INTEGER :: year,vd
	471	REAL(wp) :: s,h,p,en,p1
	472	!!
	473	INTEGER :: yr,ilc,icent,it,iday,ild,ipos,nn,iyd
[911]	474	REAL(wp) :: cycle,t,td,delt(84),delta,deltat
[1125]	475	!!
	476	DATA delt /-5.04, -3.90, -2.87, -0.58, 0.71, 1.80, &
	477	& 3.08, 4.63, 5.86, 7.21, 8.58, 10.50, 12.10, &
	478	& 12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39, &
	479	& 19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68, &
	480	& 22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63, &
	481	& 23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80, &
	482	& 24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89, &
	483	& 27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96, &
	484	& 31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39, &
	485	& 33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87, &
	486	& 38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00, &
	487	& 45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81, &
	488	& 52.57 /
	489	!!----------------------------------------------------------------------
[911]	490
	491	cycle = 360.0
	492	ilc = 0
[1125]	493	icent = year / 100
	494	yr = year - icent * 100
[911]	495	t = icent - 20
	496	!
	497	! for the following equations
	498	! time origin is fixed at 00 hr of jan 1st,2000.
	499	! see notes by cartwright
	500	!
[1125]	501	!!gm old coding style, use CASE instead and avoid GOTO (obsolescence in fortran 90)
	502	!!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90
[911]	503	it = icent - 20
	504	if (it) 1,2,2
	505	1 iday = it/4 -it
	506	go to 3
	507	2 iday = (it+3)/4 - it
	508	!
	509	! t is in julian century
	510	! correction in gegorian calander where only century year divisible
	511	! by 4 is leap year.
	512	!
	513	3 continue
	514	!
	515	td = 0.0
	516	!
[1125]	517	!!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90
[911]	518	if (yr) 4,5,4
	519	!
	520	4 iyd = 365*yr
	521	ild = (yr-1)/4
	522	if((icent - (icent/4)*4) .eq. 0) ilc = 1
	523	td = iyd + ild + ilc
	524	!
	525	5 td = td + iday + vd -1.0 - 0.5
	526	t = t + (td/36525.0)
	527	!
	528	ipos=year-1899
	529	if (ipos .lt. 0) go to 7
	530	if (ipos .gt. 83) go to 6
	531	!
	532	delta = (delt(ipos+1)+delt(ipos))/2.0
	533	go to 7
	534	!
	535	6 delta= (65.0-50.5)/20.0*(year-1980)+50.5
	536	!
	537	7 deltat = delta * 1.0e-6
	538	!
[1125]	539	!!gm precision of the computation : example for s it should be replace by:
	540	!!gm s = 218.3165 + (481267.8813 - 0.0016t)t + 152.0*deltat ==> more precise modify the last digits results
[911]	541	s = 218.3165 + 481267.8813t - 0.0016tt + 152.0deltat
[1125]	542	h = 280.4661 + 36000.7698 t + 0.0003tt + 11.0deltat
	543	p = 83.3535 + 4069.0139 t - 0.0103t*t + deltat
	544	en = 234.9555 + 1934.1363 t - 0.0021t*t + deltat
	545	p1 = 282.9384 + 1.7195 t + 0.0005t*t
	546	!
	547	nn = s / cycle
	548	s = s - nn * cycle
	549	IF( s < 0.e0 ) s = s + cycle
	550	!
	551	nn = h / cycle
	552	h = h - cycle * nn
	553	IF( h < 0.e0 ) h = h + cycle
	554	!
	555	nn = p / cycle
	556	p = p - cycle * nn
	557	IF( p < 0.e0) p = p + cycle
	558	!
	559	nn = en / cycle
	560	en = en - cycle * nn
	561	IF( en < 0.e0 ) en = en + cycle
[911]	562	en = cycle - en
[1125]	563	!
	564	nn = p1 / cycle
	565	p1 = p1 - nn * cycle
	566	!
	567	END SUBROUTINE shpen
[911]	568
	569
[1125]	570	SUBROUTINE ufset( p, cn, b, a )
	571	!!----------------------------------------------------------------------
	572	!! * SUBROUTINE ufset *
	573	!!
[911]	574	!! ** Purpose : - calculate nodal parameters for the tides
	575	!!
[1125]	576	!!----------------------------------------------------------------------
	577	!!gm doctor naming of dummy argument variables!!! and all local variables
	578	!!gm nc is jptides_max ????
[911]	579	integer nc
	580	parameter (nc=15)
[1125]	581	REAL(wp) p,cn
	582	!!
	583
	584	!!gm rad is already a public variable defined in phycst.F90 .... ==> doctor norme local real start with "z"
	585	REAL(wp) :: w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad
	586	REAL(wp) :: a(nc), b(nc)
	587	INTEGER :: ndc, k
	588	!!----------------------------------------------------------------------
[911]	589
[1125]	590	ndc = nc
	591
[911]	592	! a=f , b =u
	593	! t is zero as compared to tifa.
[1125]	594	!! use rad defined in phycst (i.e. add a USE phycst at the begining of the module
[911]	595	rad = 6.2831852d0/360.0
[1125]	596	pw = p * rad
	597	nw = cn * rad
	598	w1 = cos( nw )
	599	w2 = cos( 2*nw )
	600	w3 = cos( 3*nw )
	601	w4 = sin( nw )
	602	w5 = sin( 2*nw )
	603	w6 = sin( 3*nw )
	604	w7 = 1. - 0.2505 * COS( 2pw ) - 0.1102 COS(2*pw-nw ) &
	605	& - 0.156 * COS( 2pw-2nw ) - 0.037 * COS( nw )
	606	w8 = - 0.2505 * SIN( 2pw ) - 0.1102 SIN(2*pw-nw ) &
	607	& - 0.0156 * SIN( 2pw-2nw ) - 0.037 * SIN( nw )
	608	!
	609	a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3
	610	b(1) = 0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6
	611	! q1
[911]	612	a(2) = a(1)
	613	b(2) = b(1)
[1125]	614	! o1
[911]	615	a(3) = 1.0
	616	b(3) = 0.0
[1125]	617	! p1
[911]	618	a(4) = 1.0
	619	b(4) = 0.0
[1125]	620	! s1
[911]	621	a(5) = 1.0060+0.1150w1- 0.0088w2 +0.0006*w3
	622	b(5) = -0.1546w4 + 0.0119w5 -0.0012*w6
[1125]	623	! k1
[911]	624	a(6) =1.0004 -0.0373w1+ 0.0002w2
	625	b(6) = -0.0374*w4
[1125]	626	! 2n2
[911]	627	a(7) = a(6)
	628	b(7) = b(6)
[1125]	629	! mu2
[911]	630	a(8) = a(6)
	631	b(8) = b(6)
[1125]	632	! n2
[911]	633	a(9) = a(6)
	634	b(9) = b(6)
[1125]	635	! nu2
[911]	636	a(10) = a(6)
	637	b(10) = b(6)
[1125]	638	! m2
	639	a(11) = SQRT( w7 * w7 + w8 * w8 )
	640	b(11) = ATAN( w8 / w7 )
	641	!!gmuse rpi instead of 3.141992 ??? true pi is rpi=3.141592653589793_wp ..... ????
	642	IF( w7 < 0.e0 ) b(11) = b(11) + 3.141992
	643	! l2
[911]	644	a(12) = 1.0
	645	b(12) = 0.0
[1125]	646	! t2
[911]	647	a(13)= a(12)
	648	b(13)= b(12)
[1125]	649	! s2
[911]	650	a(14) = 1.0241+0.2863w1+0.0083w2 -0.0015*w3
	651	b(14) = -0.3096w4 + 0.0119w5 - 0.0007*w6
[1125]	652	! k2
[911]	653	a(15) = a(6)*a(6)
	654	b(15) = 2*b(6)
[1125]	655	! m4
	656	!!gm old coding, remove GOTO and label of lines
	657	!!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90
	658	DO 40 k = 1,ndc
	659	b(k) = b(k)/rad
	660	32 if (b(k)) 34,35,35
	661	34 b(k) = b(k) + 360.0
	662	go to 32
	663	35 if (b(k)-360.0) 40,37,37
	664	37 b(k) = b(k)-360.0
	665	go to 35
[911]	666	40 continue
[1125]	667	!
	668	END SUBROUTINE ufset
	669
[911]	670
[1125]	671	SUBROUTINE vset( s,h,p,en,p1,v)
	672	!!----------------------------------------------------------------------
	673	!! * SUBROUTINE vset *
	674	!!
[911]	675	!! ** Purpose : - calculate tidal phases for 0000gmt on start day of run
	676	!!
[1125]	677	!!----------------------------------------------------------------------
	678	!!gm doctor naming of dummy argument variables!!! and all local variables
	679	!!gm nc is jptides_max ????
	680	!!gm en argument is not used: suppress it ?
[911]	681	integer nc
	682	parameter (nc=15)
	683	real(wp) s,h,p,en,p1
	684	real(wp) v(nc)
[1125]	685	!!
	686	integer ndc, k
	687	!!----------------------------------------------------------------------
[911]	688
	689	ndc = nc
[1125]	690	! v s are computed here.
	691	v(1) =-3*s +h +p +270 ! Q1
	692	v(2) =-2*s +h +270.0 ! O1
	693	v(3) =-h +270 ! P1
	694	v(4) =180 ! S1
	695	v(5) =h +90.0 ! K1
	696	v(6) =-4s +2h +2*p ! 2N2
	697	v(7) =-4*(s-h) ! MU2
	698	v(8) =-3s +2h +p ! N2
	699	v(9) =-3s +4h -p ! MU2
	700	v(10) =-2s +2h ! M2
	701	v(11) =-s +2*h -p +180 ! L2
	702	v(12) =-h +p1 ! T2
	703	v(13) =0 ! S2
	704	v(14) =h+h ! K2
	705	v(15) =2*v(10) ! M4
[911]	706	!
[1125]	707	!!gm old coding, remove GOTO and label of lines
	708	!!gm obsol( 1): Arithmetic IF statement is used ===> remove this in Fortran 90
[911]	709	do 72 k = 1, ndc
[1125]	710	69 if( v(k) ) 70,71,71
	711	70 v(k) = v(k)+360.0
	712	go to 69
	713	71 if( v(k) - 360.0 ) 72,73,73
	714	73 v(k) = v(k)-360.0
	715	go to 71
[911]	716	72 continue
[1125]	717	!
	718	END SUBROUTINE vset
[911]	719
	720	#else
[1125]	721	!!----------------------------------------------------------------------
	722	!! Dummy module NO Unstruct Open Boundary Conditions for tides
	723	!!----------------------------------------------------------------------
	724	!!gm are you sure we need to define filtide and tide_cpt ?
	725	CHARACTER(len=80), PUBLIC :: filtide !: Filename root for tidal input files
	726	CHARACTER(len=4 ), PUBLIC, DIMENSION(1) :: tide_cpt !: Names of tidal components used.
[911]	727
	728	CONTAINS
[1125]	729	SUBROUTINE tide_init ! Empty routine
[911]	730	END SUBROUTINE tide_init
[1125]	731	SUBROUTINE tide_data ! Empty routine
[911]	732	END SUBROUTINE tide_data
[1125]	733	SUBROUTINE tide_update( kt, kit ) ! Empty routine
	734	WRITE(,) 'tide_update: You should not have seen this print! error?', kt, kit
[911]	735	END SUBROUTINE tide_update
	736	#endif
	737
[1125]	738	!!======================================================================
[911]	739	END MODULE bdytides

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

Download in other formats: