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.

tide_mod.F90 in NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE/tide_mod.F90 @ 10793

Last change on this file since 10793 was 10793, checked in by smueller, 5 years ago
Replacement of a hard-coded coefficient by namelist parameter rn_tide_gamma (ticket #2194)
Property svn:keywords set to `Id`
File size: 26.7 KB

Rev	Line
[2956]	1	MODULE tide_mod
[4292]	2	!!======================================================================
	3	!! * MODULE tide_mod *
	4	!! Compute nodal modulations corrections and pulsations
	5	!!======================================================================
	6	!! History : 1.0 ! 2007 (O. Le Galloudec) Original code
	7	!!----------------------------------------------------------------------
[10772]	8	USE oce ! ocean dynamics and tracers variables
[4292]	9	USE dom_oce ! ocean space and time domain
	10	USE phycst ! physical constant
	11	USE daymod ! calendar
[10772]	12	!
	13	USE in_out_manager ! I/O units
	14	USE iom ! xIOs server
	15	USE ioipsl ! NetCDF IPSL library
	16	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[2956]	17
[4292]	18	IMPLICIT NONE
	19	PRIVATE
[2956]	20
[10772]	21	PUBLIC tide_init
[4292]	22	PUBLIC tide_harmo ! called by tideini and diaharm modules
	23	PUBLIC tide_init_Wave ! called by tideini and diaharm modules
[10773]	24	PUBLIC tide_init_load
	25	PUBLIC tide_init_potential
[10777]	26	PUBLIC upd_tide ! called in dynspg_... modules
[2956]	27
[4292]	28	INTEGER, PUBLIC, PARAMETER :: jpmax_harmo = 19 !: maximum number of harmonic
[2956]	29
[4292]	30	TYPE, PUBLIC :: tide
	31	CHARACTER(LEN=4) :: cname_tide
	32	REAL(wp) :: equitide
	33	INTEGER :: nutide
	34	INTEGER :: nt, ns, nh, np, np1, shift
	35	INTEGER :: nksi, nnu0, nnu1, nnu2, R
	36	INTEGER :: nformula
	37	END TYPE tide
[2956]	38
[4292]	39	TYPE(tide), PUBLIC, DIMENSION(jpmax_harmo) :: Wave !:
[2956]	40
[10772]	41	REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: omega_tide !:
	42	REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: v0tide !:
	43	REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: utide !:
	44	REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: ftide !:
	45
	46	LOGICAL , PUBLIC :: ln_tide !:
	47	LOGICAL , PUBLIC :: ln_tide_pot !:
	48	LOGICAL , PUBLIC :: ln_read_load !:
	49	LOGICAL , PUBLIC :: ln_scal_load !:
	50	LOGICAL , PUBLIC :: ln_tide_ramp !:
	51	INTEGER , PUBLIC :: nb_harmo !:
	52	INTEGER , PUBLIC :: kt_tide !:
	53	REAL(wp), PUBLIC :: rdttideramp !:
	54	REAL(wp), PUBLIC :: rn_scal_load !:
	55	CHARACTER(lc), PUBLIC :: cn_tide_load !:
[10793]	56	REAL(wp) :: rn_tide_gamma ! Tidal tilt factor
[10772]	57
	58	INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:) :: ntide !:
	59
[10773]	60	REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: pot_astro !: tidal potential
	61	REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: amp_pot, phi_pot
	62	REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) :: amp_load, phi_load
	63
	64
[4292]	65	REAL(wp) :: sh_T, sh_s, sh_h, sh_p, sh_p1 ! astronomic angles
	66	REAL(wp) :: sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R !
	67	REAL(wp) :: sh_I, sh_x1ra, sh_N !
[2956]	68
[4292]	69	!!----------------------------------------------------------------------
[10068]	70	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5215]	71	!! $Id$
[10068]	72	!! Software governed by the CeCILL license (see ./LICENSE)
[4292]	73	!!----------------------------------------------------------------------
[2956]	74	CONTAINS
	75
[10772]	76	SUBROUTINE tide_init
	77	!!----------------------------------------------------------------------
	78	!! * ROUTINE tide_init *
	79	!!----------------------------------------------------------------------
	80	INTEGER :: ji, jk
	81	CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname
	82	INTEGER :: ios ! Local integer output status for namelist read
	83	!
[10793]	84	NAMELIST/nam_tide/ln_tide, ln_tide_pot, rn_tide_gamma, ln_scal_load, ln_read_load, cn_tide_load, &
[10772]	85	& ln_tide_ramp, rn_scal_load, rdttideramp, clname
	86	!!----------------------------------------------------------------------
	87	!
	88	! Read Namelist nam_tide
	89	REWIND( numnam_ref ) ! Namelist nam_tide in reference namelist : Tides
	90	READ ( numnam_ref, nam_tide, IOSTAT = ios, ERR = 901)
	91	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_tide in reference namelist', lwp )
	92	!
	93	REWIND( numnam_cfg ) ! Namelist nam_tide in configuration namelist : Tides
	94	READ ( numnam_cfg, nam_tide, IOSTAT = ios, ERR = 902 )
	95	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_tide in configuration namelist', lwp )
	96	IF(lwm) WRITE ( numond, nam_tide )
	97	!
	98	IF( ln_tide ) THEN
	99	IF (lwp) THEN
	100	WRITE(numout,*)
	101	WRITE(numout,*) 'tide_init : Initialization of the tidal components'
	102	WRITE(numout,*) '~~~~~~~~~ '
	103	WRITE(numout,*) ' Namelist nam_tide'
	104	WRITE(numout,*) ' Use tidal components ln_tide = ', ln_tide
	105	WRITE(numout,*) ' Apply astronomical potential ln_tide_pot = ', ln_tide_pot
[10793]	106	WRITE(numout,*) ' Tidal tilt factor rn_tide_gamma= ', rn_tide_gamma
[10772]	107	WRITE(numout,*) ' Use scalar approx. for load potential ln_scal_load = ', ln_scal_load
	108	WRITE(numout,*) ' Read load potential from file ln_read_load = ', ln_read_load
	109	WRITE(numout,*) ' Apply ramp on tides at startup ln_tide_ramp = ', ln_tide_ramp
	110	WRITE(numout,*) ' Fraction of SSH used in scal. approx. rn_scal_load = ', rn_scal_load
	111	WRITE(numout,*) ' Duration (days) of ramp rdttideramp = ', rdttideramp
	112	ENDIF
	113	ELSE
	114	rn_scal_load = 0._wp
	115
	116	IF(lwp) WRITE(numout,*)
	117	IF(lwp) WRITE(numout,*) 'tide_init : tidal components not used (ln_tide = F)'
	118	IF(lwp) WRITE(numout,*) '~~~~~~~~~ '
	119	RETURN
	120	ENDIF
	121	!
	122	CALL tide_init_Wave
	123	!
	124	nb_harmo=0
	125	DO jk = 1, jpmax_harmo
	126	DO ji = 1,jpmax_harmo
	127	IF( TRIM(clname(jk)) == Wave(ji)%cname_tide ) nb_harmo = nb_harmo + 1
	128	END DO
	129	END DO
	130	!
	131	! Ensure that tidal components have been set in namelist_cfg
	132	IF( nb_harmo == 0 ) CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' )
	133	!
	134	IF( ln_read_load.AND.(.NOT.ln_tide_pot) ) &
	135	& CALL ctl_stop('ln_read_load requires ln_tide_pot')
	136	IF( ln_scal_load.AND.(.NOT.ln_tide_pot) ) &
	137	& CALL ctl_stop('ln_scal_load requires ln_tide_pot')
	138	IF( ln_scal_load.AND.ln_read_load ) &
	139	& CALL ctl_stop('Choose between ln_scal_load and ln_read_load')
	140	IF( ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rdttideramp) ) &
	141	& CALL ctl_stop('rdttideramp must be lower than run duration')
	142	IF( ln_tide_ramp.AND.(rdttideramp<0.) ) &
	143	& CALL ctl_stop('rdttideramp must be positive')
	144	!
	145	ALLOCATE( ntide(nb_harmo) )
	146	DO jk = 1, nb_harmo
	147	DO ji = 1, jpmax_harmo
	148	IF( TRIM(clname(jk)) == Wave(ji)%cname_tide ) THEN
	149	ntide(jk) = ji
	150	EXIT
	151	ENDIF
	152	END DO
	153	END DO
	154	!
	155	ALLOCATE( omega_tide(nb_harmo), v0tide (nb_harmo), &
	156	& utide (nb_harmo), ftide (nb_harmo) )
	157	kt_tide = nit000
	158	!
	159	IF (.NOT.ln_scal_load ) rn_scal_load = 0._wp
	160	!
[10773]	161	ALLOCATE( amp_pot(jpi,jpj,nb_harmo), &
	162	& phi_pot(jpi,jpj,nb_harmo), pot_astro(jpi,jpj) )
	163	IF( ln_read_load ) THEN
	164	ALLOCATE( amp_load(jpi,jpj,nb_harmo), phi_load(jpi,jpj,nb_harmo) )
	165	ENDIF
	166	!
[10772]	167	END SUBROUTINE tide_init
	168
	169
[4292]	170	SUBROUTINE tide_init_Wave
	171	# include "tide.h90"
	172	END SUBROUTINE tide_init_Wave
[2956]	173
	174
[10773]	175	SUBROUTINE tide_init_potential
	176	!!----------------------------------------------------------------------
	177	!! * ROUTINE tide_init_potential *
	178	!!----------------------------------------------------------------------
	179	INTEGER :: ji, jj, jk ! dummy loop indices
	180	REAL(wp) :: zcons, ztmp1, ztmp2, zlat, zlon, ztmp, zamp, zcs ! local scalar
	181	!!----------------------------------------------------------------------
	182
	183	DO jk = 1, nb_harmo
[10793]	184	zcons = rn_tide_gamma * Wave(ntide(jk))%equitide * ftide(jk)
[10773]	185	DO ji = 1, jpi
	186	DO jj = 1, jpj
	187	ztmp1 = ftide(jk) * amp_pot(ji,jj,jk) * COS( phi_pot(ji,jj,jk) + v0tide(jk) + utide(jk) )
	188	ztmp2 = -ftide(jk) * amp_pot(ji,jj,jk) * SIN( phi_pot(ji,jj,jk) + v0tide(jk) + utide(jk) )
	189	zlat = gphit(ji,jj)*rad !! latitude en radian
	190	zlon = glamt(ji,jj)*rad !! longitude en radian
	191	ztmp = v0tide(jk) + utide(jk) + Wave(ntide(jk))%nutide * zlon
	192	! le potentiel est composé des effets des astres:
	193	IF ( Wave(ntide(jk))%nutide == 1 ) THEN ; zcs = zcons * SIN( 2._wp*zlat )
	194	ELSEIF( Wave(ntide(jk))%nutide == 2 ) THEN ; zcs = zcons * COS( zlat )**2
	195	ELSE ; zcs = 0._wp
	196	ENDIF
	197	ztmp1 = ztmp1 + zcs * COS( ztmp )
	198	ztmp2 = ztmp2 - zcs * SIN( ztmp )
	199	zamp = SQRT( ztmp1ztmp1 + ztmp2ztmp2 )
	200	amp_pot(ji,jj,jk) = zamp
	201	phi_pot(ji,jj,jk) = ATAN2( -ztmp2 / MAX( 1.e-10_wp , zamp ) , &
	202	& ztmp1 / MAX( 1.e-10_wp, zamp ) )
	203	END DO
	204	END DO
	205	END DO
	206	!
	207	END SUBROUTINE tide_init_potential
	208
	209
	210	SUBROUTINE tide_init_load
	211	!!----------------------------------------------------------------------
	212	!! * ROUTINE tide_init_load *
	213	!!----------------------------------------------------------------------
	214	INTEGER :: inum ! Logical unit of input file
	215	INTEGER :: ji, jj, itide ! dummy loop indices
	216	REAL(wp), DIMENSION(jpi,jpj) :: ztr, zti !: workspace to read in tidal harmonics data
	217	!!----------------------------------------------------------------------
	218	IF(lwp) THEN
	219	WRITE(numout,*)
	220	WRITE(numout,*) 'tide_init_load : Initialization of load potential from file'
	221	WRITE(numout,*) '~~~~~~~~~~~~~~ '
	222	ENDIF
	223	!
	224	CALL iom_open ( cn_tide_load , inum )
	225	!
	226	DO itide = 1, nb_harmo
	227	CALL iom_get ( inum, jpdom_data,TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
	228	CALL iom_get ( inum, jpdom_data,TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) )
	229	!
	230	DO ji=1,jpi
	231	DO jj=1,jpj
	232	amp_load(ji,jj,itide) = SQRT( ztr(ji,jj)2. + zti(ji,jj)2. )
	233	phi_load(ji,jj,itide) = ATAN2(-zti(ji,jj), ztr(ji,jj) )
	234	END DO
	235	END DO
	236	!
	237	END DO
	238	CALL iom_close( inum )
	239	!
	240	END SUBROUTINE tide_init_load
	241
	242
[4292]	243	SUBROUTINE tide_harmo( pomega, pvt, put , pcor, ktide ,kc)
	244	!!----------------------------------------------------------------------
	245	!!----------------------------------------------------------------------
	246	INTEGER , DIMENSION(kc), INTENT(in ) :: ktide ! Indice of tidal constituents
	247	INTEGER , INTENT(in ) :: kc ! Total number of tidal constituents
	248	REAL(wp), DIMENSION(kc), INTENT(out) :: pomega ! pulsation in radians/s
	249	REAL(wp), DIMENSION(kc), INTENT(out) :: pvt, put, pcor !
	250	!!----------------------------------------------------------------------
	251	!
	252	CALL astronomic_angle
	253	CALL tide_pulse( pomega, ktide ,kc )
	254	CALL tide_vuf ( pvt, put, pcor, ktide ,kc )
	255	!
	256	END SUBROUTINE tide_harmo
[2956]	257
	258
[4292]	259	SUBROUTINE astronomic_angle
	260	!!----------------------------------------------------------------------
	261	!! tj is time elapsed since 1st January 1900, 0 hour, counted in julian
	262	!! century (e.g. time in days divide by 36525)
	263	!!----------------------------------------------------------------------
	264	REAL(wp) :: cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2
	265	REAL(wp) :: zqy , zsy, zday, zdj, zhfrac
	266	!!----------------------------------------------------------------------
	267	!
	268	zqy = AINT( (nyear-1901.)/4. )
	269	zsy = nyear - 1900.
	270	!
	271	zdj = dayjul( nyear, nmonth, nday )
	272	zday = zdj + zqy - 1.
	273	!
	274	zhfrac = nsec_day / 3600.
	275	!
	276	!----------------------------------------------------------------------
	277	! Sh_n Longitude of ascending lunar node
	278	!----------------------------------------------------------------------
	279	sh_N=(259.1560564-19.328185764zsy-.0529539336zday-.0022064139zhfrac)rad
	280	!----------------------------------------------------------------------
	281	! T mean solar angle (Greenwhich time)
	282	!----------------------------------------------------------------------
	283	sh_T=(180.+zhfrac(360./24.))rad
	284	!----------------------------------------------------------------------
	285	! h mean solar Longitude
	286	!----------------------------------------------------------------------
	287	sh_h=(280.1895014-.238724988zsy+.9856473288zday+.0410686387zhfrac)rad
	288	!----------------------------------------------------------------------
	289	! s mean lunar Longitude
	290	!----------------------------------------------------------------------
	291	sh_s=(277.0256206+129.38482032zsy+13.176396768zday+.549016532zhfrac)rad
	292	!----------------------------------------------------------------------
	293	! p1 Longitude of solar perigee
	294	!----------------------------------------------------------------------
	295	sh_p1=(281.2208569+.01717836zsy+.000047064zday+.000001961zhfrac)rad
	296	!----------------------------------------------------------------------
	297	! p Longitude of lunar perigee
	298	!----------------------------------------------------------------------
	299	sh_p=(334.3837214+40.66246584zsy+.111404016zday+.004641834zhfrac)rad
[2956]	300
[4292]	301	sh_N = MOD( sh_N ,2*rpi )
	302	sh_s = MOD( sh_s ,2*rpi )
	303	sh_h = MOD( sh_h, 2*rpi )
	304	sh_p = MOD( sh_p, 2*rpi )
	305	sh_p1= MOD( sh_p1,2*rpi )
[2956]	306
[4292]	307	cosI = 0.913694997 -0.035692561 *cos(sh_N)
[2956]	308
[4292]	309	sh_I = ACOS( cosI )
[2956]	310
[4292]	311	sin2I = sin(sh_I)
	312	sh_tgn2 = tan(sh_N/2.0)
[2956]	313
[4292]	314	at1=atan(1.01883*sh_tgn2)
	315	at2=atan(0.64412*sh_tgn2)
[2956]	316
[4292]	317	sh_xi=-at1-at2+sh_N
[2956]	318
[4292]	319	IF( sh_N > rpi ) sh_xi=sh_xi-2.0*rpi
[2956]	320
[4292]	321	sh_nu = at1 - at2
[2956]	322
[4292]	323	!----------------------------------------------------------------------
	324	! For constituents l2 k1 k2
	325	!----------------------------------------------------------------------
[2956]	326
[4292]	327	tgI2 = tan(sh_I/2.0)
	328	P1 = sh_p-sh_xi
[2956]	329
[4292]	330	t2 = tgI2*tgI2
	331	t4 = t2*t2
	332	sh_x1ra = sqrt( 1.0-12.0t2cos(2.0P1)+36.0t4 )
[2956]	333
[4292]	334	p = sin(2.0*P1)
	335	q = 1.0/(6.0t2)-cos(2.0P1)
	336	sh_R = atan(p/q)
[2956]	337
[4292]	338	p = sin(2.0sh_I)sin(sh_nu)
	339	q = sin(2.0sh_I)cos(sh_nu)+0.3347
	340	sh_nuprim = atan(p/q)
[2956]	341
[4292]	342	s2 = sin(sh_I)*sin(sh_I)
	343	p = s2sin(2.0sh_nu)
	344	q = s2cos(2.0sh_nu)+0.0727
	345	sh_nusec = 0.5*atan(p/q)
	346	!
	347	END SUBROUTINE astronomic_angle
[2956]	348
	349
[4292]	350	SUBROUTINE tide_pulse( pomega, ktide ,kc )
	351	!!----------------------------------------------------------------------
	352	!! * ROUTINE tide_pulse *
	353	!!
	354	!! ** Purpose : Compute tidal frequencies
	355	!!----------------------------------------------------------------------
	356	INTEGER , INTENT(in ) :: kc ! Total number of tidal constituents
	357	INTEGER , DIMENSION(kc), INTENT(in ) :: ktide ! Indice of tidal constituents
	358	REAL(wp), DIMENSION(kc), INTENT(out) :: pomega ! pulsation in radians/s
	359	!
	360	INTEGER :: jh
	361	REAL(wp) :: zscale
	362	REAL(wp) :: zomega_T = 13149000.0_wp
	363	REAL(wp) :: zomega_s = 481267.892_wp
	364	REAL(wp) :: zomega_h = 36000.76892_wp
	365	REAL(wp) :: zomega_p = 4069.0322056_wp
	366	REAL(wp) :: zomega_n = 1934.1423972_wp
	367	REAL(wp) :: zomega_p1= 1.719175_wp
	368	!!----------------------------------------------------------------------
	369	!
	370	zscale = rad / ( 36525._wp * 86400._wp )
	371	!
	372	DO jh = 1, kc
	373	pomega(jh) = ( zomega_T * Wave( ktide(jh) )%nT &
	374	& + zomega_s * Wave( ktide(jh) )%ns &
	375	& + zomega_h * Wave( ktide(jh) )%nh &
	376	& + zomega_p * Wave( ktide(jh) )%np &
	377	& + zomega_p1* Wave( ktide(jh) )%np1 ) * zscale
	378	END DO
	379	!
	380	END SUBROUTINE tide_pulse
[2956]	381
	382
[4292]	383	SUBROUTINE tide_vuf( pvt, put, pcor, ktide ,kc )
	384	!!----------------------------------------------------------------------
	385	!! * ROUTINE tide_vuf *
	386	!!
	387	!! ** Purpose : Compute nodal modulation corrections
	388	!!
	389	!! ** Outputs : vt: Phase of tidal potential relative to Greenwich (radians)
	390	!! ut: Phase correction u due to nodal motion (radians)
	391	!! ft: Nodal correction factor
	392	!!----------------------------------------------------------------------
	393	INTEGER , INTENT(in ) :: kc ! Total number of tidal constituents
	394	INTEGER , DIMENSION(kc), INTENT(in ) :: ktide ! Indice of tidal constituents
	395	REAL(wp), DIMENSION(kc), INTENT(out) :: pvt, put, pcor !
	396	!
	397	INTEGER :: jh ! dummy loop index
	398	!!----------------------------------------------------------------------
	399	!
	400	DO jh = 1, kc
	401	! Phase of the tidal potential relative to the Greenwhich
	402	! meridian (e.g. the position of the fictuous celestial body). Units are radian:
	403	pvt(jh) = sh_T * Wave( ktide(jh) )%nT &
	404	& + sh_s * Wave( ktide(jh) )%ns &
	405	& + sh_h * Wave( ktide(jh) )%nh &
	406	& + sh_p * Wave( ktide(jh) )%np &
	407	& + sh_p1* Wave( ktide(jh) )%np1 &
	408	& + Wave( ktide(jh) )%shift * rad
	409	!
	410	! Phase correction u due to nodal motion. Units are radian:
	411	put(jh) = sh_xi * Wave( ktide(jh) )%nksi &
	412	& + sh_nu * Wave( ktide(jh) )%nnu0 &
	413	& + sh_nuprim * Wave( ktide(jh) )%nnu1 &
	414	& + sh_nusec * Wave( ktide(jh) )%nnu2 &
	415	& + sh_R * Wave( ktide(jh) )%R
[2956]	416
[4292]	417	! Nodal correction factor:
	418	pcor(jh) = nodal_factort( Wave( ktide(jh) )%nformula )
	419	END DO
	420	!
	421	END SUBROUTINE tide_vuf
[2956]	422
	423
[4292]	424	RECURSIVE FUNCTION nodal_factort( kformula ) RESULT( zf )
	425	!!----------------------------------------------------------------------
	426	!!----------------------------------------------------------------------
	427	INTEGER, INTENT(in) :: kformula
	428	!
	429	REAL(wp) :: zf
	430	REAL(wp) :: zs, zf1, zf2
	431	!!----------------------------------------------------------------------
	432	!
	433	SELECT CASE( kformula )
	434	!
	435	CASE( 0 ) !== formule 0, solar waves
	436	zf = 1.0
	437	!
	438	CASE( 1 ) !== formule 1, compound waves (78 x 78)
	439	zf=nodal_factort(78)
	440	zf = zf * zf
	441	!
	442	CASE ( 2 ) !== formule 2, compound waves (78 x 0) === (78)
	443	zf1= nodal_factort(78)
	444	zf = nodal_factort( 0)
	445	zf = zf1 * zf
[2956]	446	!
[4292]	447	CASE ( 4 ) !== formule 4, compound waves (78 x 235)
	448	zf1 = nodal_factort( 78)
	449	zf = nodal_factort(235)
	450	zf = zf1 * zf
	451	!
	452	CASE ( 5 ) !== formule 5, compound waves (78 *78 x 235)
	453	zf1 = nodal_factort( 78)
	454	zf = nodal_factort(235)
	455	zf = zf * zf1 * zf1
	456	!
	457	CASE ( 6 ) !== formule 6, compound waves (78 *78 x 0)
	458	zf1 = nodal_factort(78)
	459	zf = nodal_factort( 0)
	460	zf = zf * zf1 * zf1
	461	!
	462	CASE( 7 ) !== formule 7, compound waves (75 x 75)
	463	zf = nodal_factort(75)
	464	zf = zf * zf
	465	!
	466	CASE( 8 ) !== formule 8, compound waves (78 x 0 x 235)
	467	zf = nodal_factort( 78)
	468	zf1 = nodal_factort( 0)
	469	zf2 = nodal_factort(235)
	470	zf = zf * zf1 * zf2
	471	!
	472	CASE( 9 ) !== formule 9, compound waves (78 x 0 x 227)
	473	zf = nodal_factort( 78)
	474	zf1 = nodal_factort( 0)
	475	zf2 = nodal_factort(227)
	476	zf = zf * zf1 * zf2
	477	!
	478	CASE( 10 ) !== formule 10, compound waves (78 x 227)
	479	zf = nodal_factort( 78)
	480	zf1 = nodal_factort(227)
	481	zf = zf * zf1
	482	!
	483	CASE( 11 ) !== formule 11, compound waves (75 x 0)
	484	!!gm bug???? zf 2 fois !
	485	zf = nodal_factort(75)
[9023]	486	zf1 = nodal_factort( 0)
[4292]	487	zf = zf * zf1
	488	!
	489	CASE( 12 ) !== formule 12, compound waves (78 x 78 x 78 x 0)
	490	zf1 = nodal_factort(78)
	491	zf = nodal_factort( 0)
	492	zf = zf * zf1 * zf1 * zf1
	493	!
	494	CASE( 13 ) !== formule 13, compound waves (78 x 75)
	495	zf1 = nodal_factort(78)
	496	zf = nodal_factort(75)
	497	zf = zf * zf1
	498	!
	499	CASE( 14 ) !== formule 14, compound waves (235 x 0) === (235)
	500	zf = nodal_factort(235)
	501	zf1 = nodal_factort( 0)
	502	zf = zf * zf1
	503	!
	504	CASE( 15 ) !== formule 15, compound waves (235 x 75)
	505	zf = nodal_factort(235)
	506	zf1 = nodal_factort( 75)
	507	zf = zf * zf1
	508	!
	509	CASE( 16 ) !== formule 16, compound waves (78 x 0 x 0) === (78)
	510	zf = nodal_factort(78)
	511	zf1 = nodal_factort( 0)
	512	zf = zf * zf1 * zf1
	513	!
	514	CASE( 17 ) !== formule 17, compound waves (227 x 0)
	515	zf1 = nodal_factort(227)
	516	zf = nodal_factort( 0)
	517	zf = zf * zf1
	518	!
	519	CASE( 18 ) !== formule 18, compound waves (78 x 78 x 78 )
	520	zf1 = nodal_factort(78)
	521	zf = zf1 * zf1 * zf1
	522	!
	523	CASE( 19 ) !== formule 19, compound waves (78 x 0 x 0 x 0) === (78)
	524	!!gm bug2 ==>>> here identical to formule 16, a third multiplication by zf1 is missing
	525	zf = nodal_factort(78)
	526	zf1 = nodal_factort( 0)
	527	zf = zf * zf1 * zf1
	528	!
	529	CASE( 73 ) !== formule 73
	530	zs = sin(sh_I)
	531	zf = (2./3.-zs*zs)/0.5021
	532	!
	533	CASE( 74 ) !== formule 74
	534	zs = sin(sh_I)
	535	zf = zs * zs / 0.1578
	536	!
	537	CASE( 75 ) !== formule 75
	538	zs = cos(sh_I/2)
	539	zf = sin(sh_I) * zs * zs / 0.3800
	540	!
	541	CASE( 76 ) !== formule 76
	542	zf = sin(2*sh_I) / 0.7214
	543	!
	544	CASE( 77 ) !== formule 77
	545	zs = sin(sh_I/2)
	546	zf = sin(sh_I) * zs * zs / 0.0164
	547	!
	548	CASE( 78 ) !== formule 78
	549	zs = cos(sh_I/2)
	550	zf = zs * zs * zs * zs / 0.9154
	551	!
	552	CASE( 79 ) !== formule 79
	553	zs = sin(sh_I)
	554	zf = zs * zs / 0.1565
	555	!
	556	CASE( 144 ) !== formule 144
	557	zs = sin(sh_I/2)
	558	zf = ( 1-10zszs+15zszszszs ) * cos(sh_I/2) / 0.5873
	559	!
	560	CASE( 149 ) !== formule 149
	561	zs = cos(sh_I/2)
	562	zf = zszszszszs*zs / 0.8758
	563	!
	564	CASE( 215 ) !== formule 215
	565	zs = cos(sh_I/2)
	566	zf = zszszszs / 0.9154 sh_x1ra
	567	!
	568	CASE( 227 ) !== formule 227
	569	zs = sin(2*sh_I)
	570	zf = sqrt( 0.8965zszs+0.6001zscos (sh_nu)+0.1006 )
	571	!
	572	CASE ( 235 ) !== formule 235
	573	zs = sin(sh_I)
	574	zf = sqrt( 19.0444zszszszs + 2.7702zszscos(2sh_nu) + .0981 )
	575	!
	576	END SELECT
	577	!
	578	END FUNCTION nodal_factort
[2956]	579
	580
[4292]	581	FUNCTION dayjul( kyr, kmonth, kday )
	582	!!----------------------------------------------------------------------
	583	!! *** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE)
	584	!!----------------------------------------------------------------------
	585	INTEGER,INTENT(in) :: kyr, kmonth, kday
	586	!
	587	INTEGER,DIMENSION(12) :: idayt, idays
	588	INTEGER :: inc, ji
	589	REAL(wp) :: dayjul, zyq
	590	!
	591	DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./
	592	!!----------------------------------------------------------------------
	593	!
	594	idays(1) = 0.
	595	idays(2) = 31.
	596	inc = 0.
	597	zyq = MOD( kyr-1900. , 4. )
	598	IF( zyq == 0.) inc = 1.
	599	DO ji = 3, 12
	600	idays(ji)=idayt(ji)+inc
	601	END DO
	602	dayjul = idays(kmonth) + kday
	603	!
	604	END FUNCTION dayjul
[2956]	605
[10777]	606
	607	SUBROUTINE upd_tide( kt, kit, time_offset )
	608	!!----------------------------------------------------------------------
	609	!! * ROUTINE upd_tide *
	610	!!
	611	!! ** Purpose : provide at each time step the astronomical potential
	612	!!
	613	!! ** Method : computed from pulsation and amplitude of all tide components
	614	!!
	615	!! ** Action : pot_astro actronomical potential
	616	!!----------------------------------------------------------------------
	617	INTEGER, INTENT(in) :: kt ! ocean time-step index
	618	INTEGER, INTENT(in), OPTIONAL :: kit ! external mode sub-time-step index (lk_dynspg_ts=T)
	619	INTEGER, INTENT(in), OPTIONAL :: time_offset ! time offset in number
	620	! of internal steps (lk_dynspg_ts=F)
	621	! of external steps (lk_dynspg_ts=T)
	622	!
	623	INTEGER :: joffset ! local integer
	624	INTEGER :: ji, jj, jk ! dummy loop indices
	625	REAL(wp) :: zt, zramp ! local scalar
	626	REAL(wp), DIMENSION(nb_harmo) :: zwt
	627	!!----------------------------------------------------------------------
	628	!
	629	! ! tide pulsation at model time step (or sub-time-step)
	630	zt = ( kt - kt_tide ) * rdt
	631	!
	632	joffset = 0
	633	IF( PRESENT( time_offset ) ) joffset = time_offset
	634	!
	635	IF( PRESENT( kit ) ) THEN
	636	zt = zt + ( kit + joffset - 1 ) * rdt / REAL( nn_baro, wp )
	637	ELSE
	638	zt = zt + joffset * rdt
	639	ENDIF
	640	!
	641	zwt(:) = omega_tide(:) * zt
	642
	643	pot_astro(:,:) = 0._wp ! update tidal potential (sum of all harmonics)
	644	DO jk = 1, nb_harmo
	645	pot_astro(:,:) = pot_astro(:,:) + amp_pot(:,:,jk) * COS( zwt(jk) + phi_pot(:,:,jk) )
	646	END DO
	647	!
	648	IF( ln_tide_ramp ) THEN ! linear increase if asked
	649	zt = ( kt - nit000 ) * rdt
	650	IF( PRESENT( kit ) ) zt = zt + ( kit + joffset -1) * rdt / REAL( nn_baro, wp )
	651	zramp = MIN( MAX( zt / (rdttideramp*rday) , 0._wp ) , 1._wp )
	652	pot_astro(:,:) = zramp * pot_astro(:,:)
	653	ENDIF
	654	!
	655	END SUBROUTINE upd_tide
	656
[4292]	657	!!======================================================================
[2956]	658	END MODULE tide_mod

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

Download in other formats: