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

Last change on this file since 12045 was 12045, checked in by smueller, 4 years ago
Modifications related to coding style and conventions, and addition of a reference in module tide_mod (ticket #2194)
Property svn:keywords set to `Id`
File size: 34.5 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	!!----------------------------------------------------------------------
[11855]	8	USE oce, ONLY : sshn ! sea-surface height
	9	USE par_oce ! ocean parameters
	10	USE phycst, ONLY : rpi ! pi
	11	USE daymod, ONLY : ndt05 ! half-length of time step
	12	USE in_out_manager ! I/O units
	13	USE iom ! xIOs server
[2956]	14
[4292]	15	IMPLICIT NONE
	16	PRIVATE
[2956]	17
[10772]	18	PUBLIC tide_init
[10852]	19	PUBLIC tide_update ! called by stp
[10822]	20	PUBLIC tide_init_harmonics ! called internally and by module diaharm
[10777]	21	PUBLIC upd_tide ! called in dynspg_... modules
[2956]	22
[10811]	23	INTEGER, PUBLIC, PARAMETER :: jpmax_harmo = 64 !: maximum number of harmonic components
[2956]	24
[10852]	25	TYPE :: tide
[10811]	26	CHARACTER(LEN=4) :: cname_tide = ''
[4292]	27	REAL(wp) :: equitide
	28	INTEGER :: nt, ns, nh, np, np1, shift
	29	INTEGER :: nksi, nnu0, nnu1, nnu2, R
	30	INTEGER :: nformula
	31	END TYPE tide
[2956]	32
[10852]	33	TYPE(tide), DIMENSION(:), POINTER :: tide_components !: Array of selected tidal component parameters
[2956]	34
[10822]	35	TYPE, PUBLIC :: tide_harmonic !: Oscillation parameters of harmonic tidal components
	36	CHARACTER(LEN=4) :: cname_tide ! Name of component
	37	REAL(wp) :: equitide ! Amplitude of equilibrium tide
	38	REAL(wp) :: f ! Node factor
	39	REAL(wp) :: omega ! Angular velocity
	40	REAL(wp) :: v0 ! Initial phase at prime meridian
	41	REAL(wp) :: u ! Phase correction
	42	END type tide_harmonic
[10772]	43
[10822]	44	TYPE(tide_harmonic), PUBLIC, DIMENSION(:), POINTER :: tide_harmonics !: Oscillation parameters of selected tidal components
	45
[10772]	46	LOGICAL , PUBLIC :: ln_tide !:
	47	LOGICAL , PUBLIC :: ln_tide_pot !:
[11663]	48	INTEGER :: nn_tide_var ! Variant of tidal parameter set and tide-potential computation
[11764]	49	LOGICAL :: ln_tide_dia ! Enable tidal diagnostic output
[10852]	50	LOGICAL :: ln_read_load !:
[10772]	51	LOGICAL , PUBLIC :: ln_scal_load !:
	52	LOGICAL , PUBLIC :: ln_tide_ramp !:
[10811]	53	INTEGER , PUBLIC :: nb_harmo !: Number of active tidal components
[10800]	54	REAL(wp), PUBLIC :: rn_tide_ramp_dt !:
[10772]	55	REAL(wp), PUBLIC :: rn_scal_load !:
	56	CHARACTER(lc), PUBLIC :: cn_tide_load !:
[10793]	57	REAL(wp) :: rn_tide_gamma ! Tidal tilt factor
[10772]	58
[10773]	59	REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) :: pot_astro !: tidal potential
[11764]	60	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: pot_astro_comp ! tidal-potential component
[10852]	61	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: amp_pot, phi_pot
	62	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: amp_load, phi_load
[10773]	63
[10860]	64	REAL(wp) :: rn_tide_ramp_t ! Elapsed time in seconds
[10773]	65
[4292]	66	REAL(wp) :: sh_T, sh_s, sh_h, sh_p, sh_p1 ! astronomic angles
	67	REAL(wp) :: sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R !
	68	REAL(wp) :: sh_I, sh_x1ra, sh_N !
[2956]	69
[12042]	70	! Longitudes on 1 Jan 1900, 00h and angular velocities (units of deg and
	71	! deg/h, respectively. The values of these module variables have been copied
	72	! from subroutine astronomic_angle of the version of this module used in
	73	! version 4.0 of NEMO
	74	REAL(wp) :: rlon00_N = 259.1560564_wp ! Longitude of ascending lunar node
	75	REAL(wp) :: romega_N = -.0022064139_wp
	76	REAL(wp) :: rlon00_T = 180.0_wp ! Mean solar angle (GMT)
	77	REAL(wp) :: romega_T = 15.0_wp
	78	REAL(wp) :: rlon00_h = 280.1895014_wp ! Mean solar Longitude
	79	REAL(wp) :: romega_h = .0410686387_wp
	80	REAL(wp) :: rlon00_s = 277.0256206_wp ! Mean lunar Longitude
	81	REAL(wp) :: romega_s = .549016532_wp
	82	REAL(wp) :: rlon00_p1 = 281.2208569_wp ! Longitude of solar perigee
	83	REAL(wp) :: romega_p1 = .000001961_wp
	84	REAL(wp) :: rlon00_p = 334.3837214_wp ! Longitude of lunar perigee
	85	REAL(wp) :: romega_p = .004641834_wp
	86
[4292]	87	!!----------------------------------------------------------------------
[10068]	88	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5215]	89	!! $Id$
[10068]	90	!! Software governed by the CeCILL license (see ./LICENSE)
[4292]	91	!!----------------------------------------------------------------------
[2956]	92	CONTAINS
	93
[10772]	94	SUBROUTINE tide_init
	95	!!----------------------------------------------------------------------
	96	!! * ROUTINE tide_init *
	97	!!----------------------------------------------------------------------
	98	INTEGER :: ji, jk
[10811]	99	CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: sn_tide_cnames ! Names of selected tidal components
[10772]	100	INTEGER :: ios ! Local integer output status for namelist read
	101	!
[11764]	102	NAMELIST/nam_tide/ln_tide, nn_tide_var, ln_tide_dia, ln_tide_pot, rn_tide_gamma, &
[11663]	103	& ln_scal_load, ln_read_load, cn_tide_load, &
	104	& ln_tide_ramp, rn_scal_load, rn_tide_ramp_dt, &
	105	& sn_tide_cnames
[10772]	106	!!----------------------------------------------------------------------
	107	!
[10811]	108	! Initialise all array elements of sn_tide_cnames, as some of them
	109	! typically do not appear in namelist_ref or namelist_cfg
	110	sn_tide_cnames(:) = ''
[10772]	111	! Read Namelist nam_tide
	112	REWIND( numnam_ref ) ! Namelist nam_tide in reference namelist : Tides
	113	READ ( numnam_ref, nam_tide, IOSTAT = ios, ERR = 901)
	114	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_tide in reference namelist', lwp )
	115	!
	116	REWIND( numnam_cfg ) ! Namelist nam_tide in configuration namelist : Tides
	117	READ ( numnam_cfg, nam_tide, IOSTAT = ios, ERR = 902 )
	118	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_tide in configuration namelist', lwp )
	119	IF(lwm) WRITE ( numond, nam_tide )
	120	!
	121	IF( ln_tide ) THEN
	122	IF (lwp) THEN
	123	WRITE(numout,*)
	124	WRITE(numout,*) 'tide_init : Initialization of the tidal components'
	125	WRITE(numout,*) '~~~~~~~~~ '
	126	WRITE(numout,*) ' Namelist nam_tide'
[10800]	127	WRITE(numout,*) ' Use tidal components ln_tide = ', ln_tide
[11664]	128	WRITE(numout,*) ' Variant (1: default; 0: legacy option) nn_tide_var = ', nn_tide_var
[11764]	129	WRITE(numout,*) ' Tidal diagnostic output ln_tide_dia = ', ln_tide_dia
[10800]	130	WRITE(numout,*) ' Apply astronomical potential ln_tide_pot = ', ln_tide_pot
	131	WRITE(numout,*) ' Tidal tilt factor rn_tide_gamma = ', rn_tide_gamma
	132	WRITE(numout,*) ' Use scalar approx. for load potential ln_scal_load = ', ln_scal_load
	133	WRITE(numout,*) ' Read load potential from file ln_read_load = ', ln_read_load
	134	WRITE(numout,*) ' Apply ramp on tides at startup ln_tide_ramp = ', ln_tide_ramp
	135	WRITE(numout,*) ' Fraction of SSH used in scal. approx. rn_scal_load = ', rn_scal_load
	136	WRITE(numout,*) ' Duration (days) of ramp rn_tide_ramp_dt = ', rn_tide_ramp_dt
[10772]	137	ENDIF
	138	ELSE
	139	rn_scal_load = 0._wp
	140
	141	IF(lwp) WRITE(numout,*)
	142	IF(lwp) WRITE(numout,*) 'tide_init : tidal components not used (ln_tide = F)'
	143	IF(lwp) WRITE(numout,*) '~~~~~~~~~ '
	144	RETURN
	145	ENDIF
	146	!
	147	IF( ln_read_load.AND.(.NOT.ln_tide_pot) ) &
	148	& CALL ctl_stop('ln_read_load requires ln_tide_pot')
	149	IF( ln_scal_load.AND.(.NOT.ln_tide_pot) ) &
	150	& CALL ctl_stop('ln_scal_load requires ln_tide_pot')
	151	IF( ln_scal_load.AND.ln_read_load ) &
	152	& CALL ctl_stop('Choose between ln_scal_load and ln_read_load')
[10800]	153	IF( ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rn_tide_ramp_dt) ) &
	154	& CALL ctl_stop('rn_tide_ramp_dt must be lower than run duration')
	155	IF( ln_tide_ramp.AND.(rn_tide_ramp_dt<0.) ) &
	156	& CALL ctl_stop('rn_tide_ramp_dt must be positive')
[10772]	157	!
[10822]	158	! Initialise array used to store tidal oscillation parameters (frequency,
[10840]	159	! amplitude, phase); also retrieve and store array of information about
	160	! selected tidal components
	161	CALL tide_init_harmonics(sn_tide_cnames, tide_harmonics, tide_components)
[10822]	162	!
[10840]	163	! Number of active tidal components
	164	nb_harmo = size(tide_components)
	165	!
	166	! Ensure that tidal components have been set in namelist_cfg
	167	IF( nb_harmo == 0 ) CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' )
	168	!
[10772]	169	IF (.NOT.ln_scal_load ) rn_scal_load = 0._wp
	170	!
[10773]	171	ALLOCATE( amp_pot(jpi,jpj,nb_harmo), &
	172	& phi_pot(jpi,jpj,nb_harmo), pot_astro(jpi,jpj) )
[11764]	173	IF( ln_tide_dia ) ALLOCATE( pot_astro_comp(jpi,jpj) )
[10773]	174	IF( ln_read_load ) THEN
	175	ALLOCATE( amp_load(jpi,jpj,nb_harmo), phi_load(jpi,jpj,nb_harmo) )
[10855]	176	CALL tide_init_load
	177	amp_pot(:,:,:) = amp_load(:,:,:)
	178	phi_pot(:,:,:) = phi_load(:,:,:)
	179	ELSE
	180	amp_pot(:,:,:) = 0._wp
	181	phi_pot(:,:,:) = 0._wp
[10773]	182	ENDIF
	183	!
[10772]	184	END SUBROUTINE tide_init
	185
	186
[10811]	187	SUBROUTINE tide_init_components(pcnames, ptide_comp)
	188	!!----------------------------------------------------------------------
	189	!! * ROUTINE tide_init_components *
	190	!!
	191	!! Returns pointer to array of variables of type 'tide' that contain
	192	!! information about the selected tidal components
	193	!! ----------------------------------------------------------------------
[11864]	194	CHARACTER(LEN=4), DIMENSION(jpmax_harmo), INTENT(in) :: pcnames ! Names of selected components
	195	TYPE(tide), POINTER, DIMENSION(:), INTENT(out) :: ptide_comp ! Selected components
	196	INTEGER, ALLOCATABLE, DIMENSION(:) :: icomppos ! Indices of selected components
	197	INTEGER :: icomp, jk, jj, ji ! Miscellaneous integers
	198	LOGICAL :: llmatch ! Local variables used for
	199	INTEGER :: ic1, ic2 ! string comparison
	200	TYPE(tide), POINTER, DIMENSION(:) :: tide_components ! All available components
[10811]	201
	202	! Populate local array with information about all available tidal
	203	! components
	204	!
	205	! Note, here 'tide_components' locally overrides the global module
	206	! variable of the same name to enable the use of the global name in the
	207	! include file that contains the initialisation of elements of array
	208	! 'tide_components'
[11864]	209	ALLOCATE(tide_components(jpmax_harmo), icomppos(jpmax_harmo))
[10811]	210	! Initialise array of indices of the selected componenents
[11864]	211	icomppos(:) = 0
[10811]	212	! Include tidal component parameters for all available components
[11704]	213	IF (nn_tide_var < 1) THEN
	214	#define TIDE_VAR_0
[10811]	215	#include "tide.h90"
[11704]	216	#undef TIDE_VAR_0
	217	ELSE
	218	#include "tide.h90"
	219	END IF
[10811]	220	! Identify the selected components that are availble
[11864]	221	icomp = 0
[10811]	222	DO jk = 1, jpmax_harmo
	223	IF (TRIM(pcnames(jk)) /= '') THEN
[11864]	224	DO jj = 1, jpmax_harmo
	225	! Find matches between selected and available constituents
	226	! (ignore capitalisation unless legacy variant has been selected)
	227	IF (nn_tide_var < 1) THEN
	228	llmatch = (TRIM(pcnames(jk)) == TRIM(tide_components(jj)%cname_tide))
	229	ELSE
	230	llmatch = .TRUE.
	231	ji = MAX(LEN_TRIM(pcnames(jk)), LEN_TRIM(tide_components(jj)%cname_tide))
	232	DO WHILE (llmatch.AND.(ji > 0))
	233	ic1 = IACHAR(pcnames(jk)(ji:ji))
	234	IF ((ic1 >= 97).AND.(ic1 <= 122)) ic1 = ic1 - 32
	235	ic2 = IACHAR(tide_components(jj)%cname_tide(ji:ji))
	236	IF ((ic2 >= 97).AND.(ic2 <= 122)) ic2 = ic2 - 32
	237	llmatch = (ic1 == ic2)
	238	ji = ji - 1
	239	END DO
	240	END IF
	241	IF (llmatch) THEN
	242	! Count and record the match
	243	icomp = icomp + 1
	244	icomppos(icomp) = jj
	245	! Set the capitalisation of the tidal constituent identifier
	246	! as specified in the namelist
	247	tide_components(jj)%cname_tide = pcnames(jk)
	248	IF (lwp) WRITE(numout, '(10X,"Tidal component #",I2.2,36X,"= ",A4)') icomp, tide_components(jj)%cname_tide
[10811]	249	EXIT
	250	END IF
	251	END DO
[11864]	252	IF ((lwp).AND.(jj > jpmax_harmo)) WRITE(numout, '(10X,"Tidal component ",A4," is not available!")') pcnames(jk)
[10811]	253	END IF
	254	END DO
	255
	256	! Allocate and populate reduced list of components
[11864]	257	ALLOCATE(ptide_comp(icomp))
	258	DO jk = 1, icomp
	259	ptide_comp(jk) = tide_components(icomppos(jk))
[10811]	260	END DO
	261
	262	! Release local array of available components and list of selected
	263	! components
[11864]	264	DEALLOCATE(tide_components, icomppos)
[10811]	265
	266	END SUBROUTINE tide_init_components
[2956]	267
	268
[10840]	269	SUBROUTINE tide_init_harmonics(pcnames, ptide_harmo, ptide_comp)
[10822]	270	!!----------------------------------------------------------------------
	271	!! * ROUTINE tide_init_harmonics *
	272	!!
	273	!! Returns pointer to array of variables of type 'tide_harmonics' that
	274	!! contain oscillation parameters of the selected harmonic tidal
	275	!! components
	276	!! ----------------------------------------------------------------------
[10840]	277	CHARACTER(LEN=4), DIMENSION(jpmax_harmo), INTENT(in) :: pcnames ! Names of selected components
	278	TYPE(tide_harmonic), POINTER, DIMENSION(:) :: ptide_harmo ! Oscillation parameters of tidal components
	279	TYPE(tide), POINTER, DIMENSION(:), OPTIONAL :: ptide_comp ! Selected components
	280	TYPE(tide), POINTER, DIMENSION(:) :: ztcomp ! Selected components
[10822]	281
[10840]	282	! Retrieve information about selected tidal components
	283	! If requested, prepare tidal component array for returning
	284	IF (PRESENT(ptide_comp)) THEN
	285	CALL tide_init_components(pcnames, ptide_comp)
	286	ztcomp => ptide_comp
	287	ELSE
	288	CALL tide_init_components(pcnames, ztcomp)
	289	END IF
	290
[10822]	291	! Allocate and populate array of oscillation parameters
[10840]	292	ALLOCATE(ptide_harmo(size(ztcomp)))
	293	ptide_harmo(:)%cname_tide = ztcomp(:)%cname_tide
	294	ptide_harmo(:)%equitide = ztcomp(:)%equitide
	295	CALL tide_harmo(ztcomp, ptide_harmo)
[10822]	296
	297	END SUBROUTINE tide_init_harmonics
	298
	299
[10773]	300	SUBROUTINE tide_init_potential
	301	!!----------------------------------------------------------------------
	302	!! * ROUTINE tide_init_potential *
[11737]	303	!!
	304	!! *** Reference:
	305	!! CT71) Cartwright, D. E. and Tayler, R. J. (1971): New computations of
	306	!! the Tide-generating Potential. Geophys. J. R. astr. Soc. 23,
	307	!! pp. 45-74. DOI: 10.1111/j.1365-246X.1971.tb01803.x
	308	!!
[10773]	309	!!----------------------------------------------------------------------
	310	INTEGER :: ji, jj, jk ! dummy loop indices
	311	REAL(wp) :: zcons, ztmp1, ztmp2, zlat, zlon, ztmp, zamp, zcs ! local scalar
	312	!!----------------------------------------------------------------------
	313
[10855]	314	IF( ln_read_load ) THEN
	315	amp_pot(:,:,:) = amp_load(:,:,:)
	316	phi_pot(:,:,:) = phi_load(:,:,:)
	317	ELSE
	318	amp_pot(:,:,:) = 0._wp
	319	phi_pot(:,:,:) = 0._wp
	320	ENDIF
[10773]	321	DO jk = 1, nb_harmo
[10822]	322	zcons = rn_tide_gamma * tide_components(jk)%equitide * tide_harmonics(jk)%f
[10773]	323	DO ji = 1, jpi
	324	DO jj = 1, jpj
[10822]	325	ztmp1 = tide_harmonics(jk)%f * amp_pot(ji,jj,jk) * COS( phi_pot(ji,jj,jk) &
	326	& + tide_harmonics(jk)%v0 + tide_harmonics(jk)%u )
	327	ztmp2 = -tide_harmonics(jk)%f * amp_pot(ji,jj,jk) * SIN( phi_pot(ji,jj,jk) &
	328	& + tide_harmonics(jk)%v0 + tide_harmonics(jk)%u )
[10773]	329	zlat = gphit(ji,jj)*rad !! latitude en radian
	330	zlon = glamt(ji,jj)*rad !! longitude en radian
[11722]	331	ztmp = tide_harmonics(jk)%v0 + tide_harmonics(jk)%u + tide_components(jk)%nt * zlon
[10773]	332	! le potentiel est composé des effets des astres:
[11722]	333	SELECT CASE( tide_components(jk)%nt )
[11664]	334	CASE( 0 ) ! long-periodic tidal constituents (included unless
	335	zcs = zcons * ( 0.5_wp - 1.5_wp * SIN( zlat )**2 ) ! compatibility with original formulation is requested)
	336	IF ( nn_tide_var < 1 ) zcs = 0.0_wp
	337	CASE( 1 ) ! diurnal tidal constituents
	338	zcs = zcons * SIN( 2.0_wp*zlat )
	339	CASE( 2 ) ! semi-diurnal tidal constituents
	340	zcs = zcons * COS( zlat )**2
[11737]	341	CASE( 3 ) ! Terdiurnal tidal constituents; the colatitude-dependent
	342	zcs = zcons * COS( zlat )**3 ! factor is sin(theta)^3 (Table 2 of CT71)
[11664]	343	CASE DEFAULT ! constituents of higher frequency are not included
	344	zcs = 0.0_wp
	345	END SELECT
[10773]	346	ztmp1 = ztmp1 + zcs * COS( ztmp )
	347	ztmp2 = ztmp2 - zcs * SIN( ztmp )
	348	zamp = SQRT( ztmp1ztmp1 + ztmp2ztmp2 )
	349	amp_pot(ji,jj,jk) = zamp
	350	phi_pot(ji,jj,jk) = ATAN2( -ztmp2 / MAX( 1.e-10_wp , zamp ) , &
	351	& ztmp1 / MAX( 1.e-10_wp, zamp ) )
	352	END DO
	353	END DO
	354	END DO
	355	!
	356	END SUBROUTINE tide_init_potential
	357
	358
	359	SUBROUTINE tide_init_load
	360	!!----------------------------------------------------------------------
	361	!! * ROUTINE tide_init_load *
	362	!!----------------------------------------------------------------------
	363	INTEGER :: inum ! Logical unit of input file
	364	INTEGER :: ji, jj, itide ! dummy loop indices
	365	REAL(wp), DIMENSION(jpi,jpj) :: ztr, zti !: workspace to read in tidal harmonics data
	366	!!----------------------------------------------------------------------
	367	IF(lwp) THEN
	368	WRITE(numout,*)
	369	WRITE(numout,*) 'tide_init_load : Initialization of load potential from file'
	370	WRITE(numout,*) '~~~~~~~~~~~~~~ '
	371	ENDIF
	372	!
	373	CALL iom_open ( cn_tide_load , inum )
	374	!
	375	DO itide = 1, nb_harmo
[10811]	376	CALL iom_get ( inum, jpdom_data,TRIM(tide_components(itide)%cname_tide)//'_z1', ztr(:,:) )
	377	CALL iom_get ( inum, jpdom_data,TRIM(tide_components(itide)%cname_tide)//'_z2', zti(:,:) )
[10773]	378	!
	379	DO ji=1,jpi
	380	DO jj=1,jpj
	381	amp_load(ji,jj,itide) = SQRT( ztr(ji,jj)2. + zti(ji,jj)2. )
	382	phi_load(ji,jj,itide) = ATAN2(-zti(ji,jj), ztr(ji,jj) )
	383	END DO
	384	END DO
	385	!
	386	END DO
	387	CALL iom_close( inum )
	388	!
	389	END SUBROUTINE tide_init_load
	390
	391
[10852]	392	SUBROUTINE tide_update( kt )
	393	!!----------------------------------------------------------------------
	394	!! * ROUTINE tide_update *
	395	!!----------------------------------------------------------------------
	396	INTEGER, INTENT( in ) :: kt ! ocean time-step
	397	INTEGER :: jk ! dummy loop index
	398	!!----------------------------------------------------------------------
	399
	400	IF( nsec_day == NINT(0.5_wp * rdt) .OR. kt == nit000 ) THEN ! start a new day
	401	!
[10856]	402	CALL tide_harmo(tide_components, tide_harmonics, ndt05) ! Update oscillation parameters of tidal components for start of current day
[10852]	403	!
	404	!
	405	IF(lwp) THEN
	406	WRITE(numout,*)
	407	WRITE(numout,*) 'tide_update : Update of the components and (re)Init. the potential at kt=', kt
	408	WRITE(numout,*) '~~~~~~~~~~~ '
	409	DO jk = 1, nb_harmo
	410	WRITE(numout,*) tide_harmonics(jk)%cname_tide, tide_harmonics(jk)%u, &
	411	& tide_harmonics(jk)%f,tide_harmonics(jk)%v0, tide_harmonics(jk)%omega
	412	END DO
	413	ENDIF
	414	!
	415	IF( ln_tide_pot ) CALL tide_init_potential
	416	!
[10860]	417	rn_tide_ramp_t = (kt - nit000)*rdt ! Elapsed time in seconds
[10852]	418	ENDIF
	419	!
	420	END SUBROUTINE tide_update
	421
	422
[10856]	423	SUBROUTINE tide_harmo( ptide_comp, ptide_harmo, psec_day )
[4292]	424	!
[10822]	425	TYPE(tide), DIMENSION(:), POINTER :: ptide_comp ! Array of selected tidal component parameters
	426	TYPE(tide_harmonic), DIMENSION(:), POINTER :: ptide_harmo ! Oscillation parameters of selected tidal components
[10856]	427	INTEGER, OPTIONAL :: psec_day ! Number of seconds since the start of the current day
[10822]	428	!
[10856]	429	IF (PRESENT(psec_day)) THEN
	430	CALL astronomic_angle(psec_day)
	431	ELSE
	432	CALL astronomic_angle(nsec_day)
	433	END IF
[10822]	434	CALL tide_pulse( ptide_comp, ptide_harmo )
	435	CALL tide_vuf( ptide_comp, ptide_harmo )
[4292]	436	!
	437	END SUBROUTINE tide_harmo
[2956]	438
	439
[10856]	440	SUBROUTINE astronomic_angle(psec_day)
[4292]	441	!!----------------------------------------------------------------------
[12042]	442	!! * ROUTINE astronomic_angle *
	443	!!
	444	!! ** Purpose : Compute astronomic angles
[4292]	445	!!----------------------------------------------------------------------
[10856]	446	INTEGER :: psec_day ! Number of seconds from midnight
[4292]	447	REAL(wp) :: cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2
[12042]	448	REAL(wp) :: zqy , zsy, zday, zdj, zhfrac, zt
[4292]	449	!!----------------------------------------------------------------------
	450	!
[12042]	451	! Computation of the time from 1 Jan 1900, 00h in years
	452	zqy = AINT( (nyear - 1901.0_wp) / 4.0_wp )
	453	zsy = nyear - 1900.0_wp
[4292]	454	!
	455	zdj = dayjul( nyear, nmonth, nday )
[12042]	456	zday = zdj + zqy - 1.0_wp
[4292]	457	!
[12042]	458	zhfrac = psec_day / 3600.0_wp
[4292]	459	!
[12042]	460	zt = zsy * 365.0_wp * 24.0_wp + zday * 24.0_wp + zhfrac
	461	!
	462	! Longitude of ascending lunar node
	463	sh_N = ( rlon00_N + romega_N * zt ) * rad
	464	sh_N = MOD( sh_N, 2*rpi )
	465	! Mean solar angle (Greenwhich time)
	466	sh_T = ( rlon00_T + romega_T * zhfrac ) * rad
	467	! Mean solar Longitude
	468	sh_h = ( rlon00_h + romega_h * zt ) * rad
	469	sh_h = MOD( sh_h, 2*rpi )
	470	! Mean lunar Longitude
	471	sh_s = ( rlon00_s + romega_s * zt ) * rad
	472	sh_s = MOD( sh_s, 2*rpi )
	473	! Longitude of solar perigee
	474	sh_p1 = ( rlon00_p1 + romega_p1 * zt ) * rad
	475	sh_p1= MOD( sh_p1, 2*rpi )
	476	! Longitude of lunar perigee
	477	sh_p = ( rlon00_p + romega_p * zt ) * rad
	478	sh_p = MOD( sh_p, 2*rpi )
[2956]	479
[12045]	480	cosI = 0.913694997_wp - 0.035692561_wp * COS( sh_N )
[2956]	481
[4292]	482	sh_I = ACOS( cosI )
[2956]	483
[4292]	484	sin2I = sin(sh_I)
	485	sh_tgn2 = tan(sh_N/2.0)
[2956]	486
[12045]	487	at1 = ATAN( 1.01883_wp * sh_tgn2 )
	488	at2 = ATAN( 0.64412_wp * sh_tgn2 )
[2956]	489
[12045]	490	sh_xi = sh_N - at1 - at2
[2956]	491
[12045]	492	IF( sh_N > rpi ) sh_xi = sh_xi - 2.0_wp * rpi
[2956]	493
[4292]	494	sh_nu = at1 - at2
[2956]	495
[12045]	496	! For computation of tidal constituents L2 K1 K2
	497	tgI2 = tan( sh_I / 2.0_wp )
	498	P1 = sh_p - sh_xi
[4292]	499	!
[12045]	500	t2 = tgI2 * tgI2
	501	t4 = t2 * t2
	502	sh_x1ra = SQRT( 1.0 - 12.0 * t2 * COS( 2.0 * P1 ) + 36.0_wp * t4 )
	503	!
	504	p = SIN( 2.0_wp * P1 )
	505	q = 1.0_wp / ( 6.0_wp * t2 ) - COS( 2.0_wp * P1 )
	506	sh_R = ATAN( p / q )
	507	!
	508	p = SIN( 2.0_wp * sh_I ) * SIN( sh_nu )
	509	q = SIN( 2.0_wp * sh_I ) * COS( sh_nu ) + 0.3347_wp
	510	sh_nuprim = ATAN( p / q )
	511	!
	512	s2 = SIN( sh_I ) * SIN( sh_I )
	513	p = s2 * SIN( 2.0_wp * sh_nu )
	514	q = s2 * COS( 2.0_wp * sh_nu ) + 0.0727_wp
	515	sh_nusec = 0.5_wp * ATAN( p / q )
	516	!
[4292]	517	END SUBROUTINE astronomic_angle
[2956]	518
	519
[10822]	520	SUBROUTINE tide_pulse( ptide_comp, ptide_harmo )
[4292]	521	!!----------------------------------------------------------------------
	522	!! * ROUTINE tide_pulse *
	523	!!
	524	!! ** Purpose : Compute tidal frequencies
	525	!!----------------------------------------------------------------------
[10822]	526	TYPE(tide), DIMENSION(:), POINTER :: ptide_comp ! Array of selected tidal component parameters
	527	TYPE(tide_harmonic), DIMENSION(:), POINTER :: ptide_harmo ! Oscillation parameters of selected tidal components
[4292]	528	!
	529	INTEGER :: jh
	530	REAL(wp) :: zscale
	531	REAL(wp) :: zomega_T = 13149000.0_wp
	532	REAL(wp) :: zomega_s = 481267.892_wp
	533	REAL(wp) :: zomega_h = 36000.76892_wp
	534	REAL(wp) :: zomega_p = 4069.0322056_wp
	535	REAL(wp) :: zomega_n = 1934.1423972_wp
	536	REAL(wp) :: zomega_p1= 1.719175_wp
	537	!!----------------------------------------------------------------------
	538	!
	539	zscale = rad / ( 36525._wp * 86400._wp )
	540	!
[10839]	541	DO jh = 1, size(ptide_harmo)
[10822]	542	ptide_harmo(jh)%omega = ( zomega_T * ptide_comp( jh )%nT &
	543	& + zomega_s * ptide_comp( jh )%ns &
	544	& + zomega_h * ptide_comp( jh )%nh &
	545	& + zomega_p * ptide_comp( jh )%np &
	546	& + zomega_p1* ptide_comp( jh )%np1 ) * zscale
[4292]	547	END DO
	548	!
	549	END SUBROUTINE tide_pulse
[2956]	550
	551
[10822]	552	SUBROUTINE tide_vuf( ptide_comp, ptide_harmo )
[4292]	553	!!----------------------------------------------------------------------
	554	!! * ROUTINE tide_vuf *
	555	!!
	556	!! ** Purpose : Compute nodal modulation corrections
	557	!!
	558	!! ** Outputs : vt: Phase of tidal potential relative to Greenwich (radians)
	559	!! ut: Phase correction u due to nodal motion (radians)
	560	!! ft: Nodal correction factor
	561	!!----------------------------------------------------------------------
[10822]	562	TYPE(tide), DIMENSION(:), POINTER :: ptide_comp ! Array of selected tidal component parameters
	563	TYPE(tide_harmonic), DIMENSION(:), POINTER :: ptide_harmo ! Oscillation parameters of selected tidal components
[4292]	564	!
	565	INTEGER :: jh ! dummy loop index
	566	!!----------------------------------------------------------------------
	567	!
[10839]	568	DO jh = 1, size(ptide_harmo)
[4292]	569	! Phase of the tidal potential relative to the Greenwhich
	570	! meridian (e.g. the position of the fictuous celestial body). Units are radian:
[10822]	571	ptide_harmo(jh)%v0 = sh_T * ptide_comp( jh )%nT &
	572	& + sh_s * ptide_comp( jh )%ns &
	573	& + sh_h * ptide_comp( jh )%nh &
	574	& + sh_p * ptide_comp( jh )%np &
	575	& + sh_p1* ptide_comp( jh )%np1 &
	576	& + ptide_comp( jh )%shift * rad
[4292]	577	!
	578	! Phase correction u due to nodal motion. Units are radian:
[10822]	579	ptide_harmo(jh)%u = sh_xi * ptide_comp( jh )%nksi &
	580	& + sh_nu * ptide_comp( jh )%nnu0 &
	581	& + sh_nuprim * ptide_comp( jh )%nnu1 &
	582	& + sh_nusec * ptide_comp( jh )%nnu2 &
	583	& + sh_R * ptide_comp( jh )%R
[2956]	584
[4292]	585	! Nodal correction factor:
[10822]	586	ptide_harmo(jh)%f = nodal_factort( ptide_comp( jh )%nformula )
[4292]	587	END DO
	588	!
	589	END SUBROUTINE tide_vuf
[2956]	590
	591
[4292]	592	RECURSIVE FUNCTION nodal_factort( kformula ) RESULT( zf )
	593	!!----------------------------------------------------------------------
[12045]	594	!! * FUNCTION nodal_factort *
	595	!!
	596	!! ** Purpose : Compute amplitude correction factors due to nodal motion
	597	!!
	598	!! ** Reference :
	599	!! S58) Schureman, P. (1958): Manual of Harmonic Analysis and
	600	!! Prediction of Tides (Revised (1940) Edition (Reprinted 1958
	601	!! with corrections). Reprinted June 2001). U.S. Department of
	602	!! Commerce, Coast and Geodetic Survey Special Publication
	603	!! No. 98. Washington DC, United States Government Printing
	604	!! Office. 317 pp. DOI: 10.25607/OBP-155.
[4292]	605	!!----------------------------------------------------------------------
[12045]	606	INTEGER, INTENT(in) :: kformula
[4292]	607	!
[12045]	608	REAL(wp) :: zf
	609	REAL(wp) :: zs, zf1, zf2
	610	CHARACTER(LEN=3) :: clformula
[4292]	611	!!----------------------------------------------------------------------
	612	!
	613	SELECT CASE( kformula )
	614	!
[12045]	615	CASE( 0 ) ! Formula 0, solar waves
[4292]	616	zf = 1.0
	617	!
[12045]	618	CASE( 1 ) ! Formula 1, compound waves (78 x 78)
	619	zf=nodal_factort( 78 )
[4292]	620	zf = zf * zf
	621	!
[12045]	622	CASE ( 4 ) ! Formula 4, compound waves (78 x 235)
	623	zf1 = nodal_factort( 78 )
[4292]	624	zf = nodal_factort(235)
	625	zf = zf1 * zf
	626	!
[12045]	627	CASE( 18 ) ! Formula 18, compound waves (78 x 78 x 78 )
	628	zf1 = nodal_factort( 78 )
[4292]	629	zf = zf1 * zf1 * zf1
	630	!
[12045]	631	CASE( 20 ) ! Formula 20, compound waves ( 78 x 78 x 78 x 78 )
	632	zf1 = nodal_factort( 78 )
[11768]	633	zf = zf1 * zf1 * zf1 * zf1
	634	!
[12045]	635	CASE( 73 ) ! Formula 73 of S58
	636	zs = SIN( sh_I )
	637	zf = ( 2.0_wp / 3.0_wp - zs * zs ) / 0.5021_wp
[4292]	638	!
[12045]	639	CASE( 74 ) ! Formula 74 of S58
	640	zs = SIN(sh_I)
	641	zf = zs * zs / 0.1578_wp
[4292]	642	!
[12045]	643	CASE( 75 ) ! Formula 75 of S58
	644	zs = COS( sh_I / 2.0_wp )
	645	zf = SIN( sh_I ) * zs * zs / 0.3800_wp
[4292]	646	!
[12045]	647	CASE( 76 ) ! Formula 76 of S58
	648	zf = SIN( 2.0_wp * sh_I ) / 0.7214_wp
[4292]	649	!
[12045]	650	CASE( 78 ) ! Formula 78 of S58
	651	zs = COS( sh_I/2 )
	652	zf = zs * zs * zs * zs / 0.9154_wp
[4292]	653	!
[12045]	654	CASE( 149 ) ! Formula 149 of S58
	655	zs = COS( sh_I/2 )
	656	zf = zs * zs * zs * zs * zs * zs / 0.8758_wp
[4292]	657	!
[12045]	658	CASE( 215 ) ! Formula 215 of S58 with typo correction (0.9154 instead of 0.9145)
	659	zs = COS( sh_I/2 )
	660	zf = zs * zs * zs * zs / 0.9154_wp * sh_x1ra
[4292]	661	!
[12045]	662	CASE( 227 ) ! Formula 227 of S58
	663	zs = SIN( 2.0_wp * sh_I )
	664	zf = SQRT( 0.8965_wp * zs * zs + 0.6001_wp * zs * COS( sh_nu ) + 0.1006_wp )
[4292]	665	!
[12045]	666	CASE ( 235 ) ! Formula 235 of S58
	667	zs = SIN( sh_I )
	668	zf = SQRT( 19.0444_wp * zs * zs * zs * zs + 2.7702_wp * zs * zs * cos( 2.0_wp * sh_nu ) + 0.0981_wp )
[4292]	669	!
[11770]	670	CASE DEFAULT
	671	WRITE( clformula, '(I3)' ) kformula
	672	CALL ctl_stop('nodal_factort: formula ' // clformula // ' is not available')
[4292]	673	END SELECT
	674	!
	675	END FUNCTION nodal_factort
[2956]	676
	677
[4292]	678	FUNCTION dayjul( kyr, kmonth, kday )
	679	!!----------------------------------------------------------------------
[12045]	680	!! * FUNCTION dayjul *
	681	!!
	682	!! Purpose : compute the Julian day
[4292]	683	!!----------------------------------------------------------------------
[12045]	684	INTEGER,INTENT(in) :: kyr, kmonth, kday
[4292]	685	!
[12045]	686	INTEGER,DIMENSION(12) :: idayt = (/ 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334 /)
	687	INTEGER,DIMENSION(12) :: idays
	688	INTEGER :: inc, ji, zyq
	689	REAL(wp) :: dayjul
[4292]	690	!!----------------------------------------------------------------------
	691	!
[12045]	692	idays(1) = 0
	693	idays(2) = 31
	694	inc = 0.0_wp
	695	zyq = MOD( kyr - 1900 , 4 )
	696	IF( zyq == 0 ) inc = 1
[4292]	697	DO ji = 3, 12
[12045]	698	idays(ji) = idayt(ji) + inc
[4292]	699	END DO
[12045]	700	dayjul = REAL( idays(kmonth) + kday, KIND=wp )
[4292]	701	!
	702	END FUNCTION dayjul
[2956]	703
[10777]	704
[10860]	705	SUBROUTINE upd_tide(pdelta)
[10777]	706	!!----------------------------------------------------------------------
	707	!! * ROUTINE upd_tide *
	708	!!
	709	!! ** Purpose : provide at each time step the astronomical potential
	710	!!
	711	!! ** Method : computed from pulsation and amplitude of all tide components
	712	!!
	713	!! ** Action : pot_astro actronomical potential
	714	!!----------------------------------------------------------------------
[10860]	715	REAL, INTENT(in) :: pdelta ! Temporal offset in seconds
	716	INTEGER :: jk ! Dummy loop index
	717	REAL(wp) :: zt, zramp ! Local scalars
	718	REAL(wp), DIMENSION(nb_harmo) :: zwt ! Temporary array
[10777]	719	!!----------------------------------------------------------------------
	720	!
[10860]	721	zwt(:) = tide_harmonics(:)%omega * pdelta
[10777]	722	!
	723	IF( ln_tide_ramp ) THEN ! linear increase if asked
[10860]	724	zt = rn_tide_ramp_t + pdelta
[10800]	725	zramp = MIN( MAX( zt / (rn_tide_ramp_dt*rday) , 0._wp ) , 1._wp )
[10777]	726	ENDIF
	727	!
[11764]	728	pot_astro(:,:) = 0._wp ! update tidal potential (sum of all harmonics)
	729	DO jk = 1, nb_harmo
	730	IF ( .NOT. ln_tide_dia ) THEN
	731	pot_astro(:,:) = pot_astro(:,:) + amp_pot(:,:,jk) * COS( zwt(jk) + phi_pot(:,:,jk) )
	732	ELSE
	733	pot_astro_comp(:,:) = amp_pot(:,:,jk) * COS( zwt(jk) + phi_pot(:,:,jk) )
	734	pot_astro(:,:) = pot_astro(:,:) + pot_astro_comp(:,:)
	735	IF ( iom_use( "tide_pot_" // TRIM( tide_harmonics(jk)%cname_tide ) ) ) THEN ! Output tidal potential (incl. load potential)
	736	IF ( ln_tide_ramp ) pot_astro_comp(:,:) = zramp * pot_astro_comp(:,:)
	737	CALL iom_put( "tide_pot_" // TRIM( tide_harmonics(jk)%cname_tide ), pot_astro_comp(:,:) )
	738	END IF
	739	END IF
	740	END DO
	741	!
	742	IF ( ln_tide_ramp ) pot_astro(:,:) = zramp * pot_astro(:,:)
	743	!
	744	IF( ln_tide_dia ) THEN ! Output total tidal potential (incl. load potential)
	745	IF ( iom_use( "tide_pot" ) ) CALL iom_put( "tide_pot", pot_astro(:,:) + rn_scal_load * sshn(:,:) )
	746	END IF
	747	!
[10777]	748	END SUBROUTINE upd_tide
	749
[4292]	750	!!======================================================================
[2956]	751	END MODULE tide_mod

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

Download in other formats: