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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 3593

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

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

Download in other formats: