Context Navigation

diaharm_fast.F90 @ 15500

Last change on this file since 15500 was 15500, checked in by hadjt, 8 months ago

diaharm_fast.F90

Found an important bug. it wasn't set up for monthly cycles, so was only reseting the "seconds of the day" count once a month. Therefore, the v0tide was introducing a jump every day, so was fitting poorly.

It is still not working properly, as the amplitudes output are too high, and the phase maps are not correct.

File size: 45.2 KB

Line
1	MODULE diaharm_fast
2	!!======================================================================
3	!! * MODULE example *
4	!! Ocean physics: On line harmonic analyser
5	!!
6	!!=====================================================================
7
8	#if defined key_diaharm_fast
9
10	!!----------------------------------------------------------------------
11	!! 'key_harm_ana' : Calculate harmonic analysis
12	!!----------------------------------------------------------------------
13	!! harm_ana :
14	!! harm_ana_init :
15	!! NB: 2017-12 : add 3D harmonic analysis of velocities
16	!! integration of Maria Luneva's development
17	!! 'key_3Ddiaharm'
18	!!----------------------------------------------------------------------
19
20	USE oce ! ocean dynamics and tracers
21	USE dom_oce ! ocean space and time domain
22	USE iom
23	USE in_out_manager ! I/O units
24	USE phycst ! physical constants
25	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
26	USE bdy_oce ! ocean open boundary conditions
27	USE bdytides ! tidal bdy forcing
28	USE daymod ! calendar
29	USE tideini
30	USE restart
31	USE ioipsl, ONLY : ju2ymds ! for calendar
32	!
33	!
34	USE timing ! preformance summary
35	USE zdf_oce
36
37	IMPLICIT NONE
38	PRIVATE
39
40	!! * Routine accessibility
41	PUBLIC dia_harm_fast ! routine called in step.F90 module
42	LOGICAL, PUBLIC, PARAMETER :: lk_diaharm_fast = .TRUE. ! to be run or not
43	LOGICAL, PUBLIC :: lk_diaharm_2D ! = .TRUE. ! to run 2d
44	LOGICAL, PUBLIC :: lk_diaharm_3D ! = .TRUE. ! to run 3d
45
46	!! * Module variables
47	INTEGER, PARAMETER :: nharm_max = jpmax_harmo ! max number of harmonics to be analysed
48	INTEGER, PARAMETER :: nhm_max = 2*nharm_max+1
49	INTEGER, PARAMETER :: nvab = 2 ! number of 3D variables
50	INTEGER :: nharm
51	INTEGER :: nhm
52	INTEGER :: & !!! toto namelist (namtoto)
53	nflag = 1 ! default value of nflag
54	REAL(wp), DIMENSION(nharm_max) :: &
55	om_tide ! tidal frequencies ( rads/sec)
56	REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:) :: &
57	bzz,c,x ! work arrays
58	REAL(wp) :: cca,ssa,zm,bt,dd_cumul
59	!
60	REAL(wp), PUBLIC :: fjulday_startharm !: Julian Day since start of harmonic analysis
61	REAL(wp), PUBLIC, ALLOCATABLE,DIMENSION(:) :: anau, anav, anaf ! nodel/phase corrections used by diaharmana
62	REAL(WP), ALLOCATABLE,SAVE,DIMENSION(:,:) :: cc,a
63	!
64	INTEGER :: nvar_2d, nvar_3d !: number of 2d and 3d variables to analyse
65	INTEGER, ALLOCATABLE,DIMENSION(:) :: m_posi_2d, m_posi_3d
66
67	! Name of variables used in the restart
68	CHARACTER( LEN = 10 ), DIMENSION(5), PARAMETER :: m_varName2d = (/'ssh','u2d','v2d','ubfr','vbfr'/)
69	CHARACTER( LEN = 10 ), DIMENSION(4), PARAMETER :: m_varName3d = (/'rho','u3d','v3d','w3d'/)
70	!
71	REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,: ) :: g_cosamp2D, g_sinamp2D, g_cumul_var2D
72	REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,:,:) :: g_cosamp3D, g_sinamp3D, g_cumul_var3D
73	!
74	REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:) :: g_out2D,h_out2D ! arrays for output
75	REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: g_out3D,h_out3D ! arrays for 3D output
76	!
77	! NAMELIST
78	LOGICAL, PUBLIC :: ln_diaharm_store !: =T Stores data for harmonic Analysis
79	LOGICAL, PUBLIC :: ln_diaharm_compute !: =T Compute harmonic Analysis
80	LOGICAL, PUBLIC :: ln_diaharm_read_restart !: =T Read restart from a previous run
81	!JT
82	LOGICAL, PUBLIC :: ln_diaharm_multiyear !: =T Read restart from a previous run
83	INTEGER, PUBLIC :: nn_diaharm_multiyear !: =T Read restart from a previous run
84	LOGICAL, PUBLIC :: ln_diaharm_update_nodal_daily !: =T update the nodes every day
85	LOGICAL, PUBLIC :: ln_diaharm_fast
86
87	!JT
88	LOGICAL, PUBLIC :: ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d
89	INTEGER :: nb_ana ! Number of harmonics to analyse
90	CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: tname ! Names of tidal constituents ('M2', 'K1',...)
91	INTEGER , ALLOCATABLE, DIMENSION(:) :: ntide_all ! INDEX within the full set of constituents (tide.h90)
92	INTEGER , ALLOCATABLE, DIMENSION(:) :: ntide_sub ! INDEX within the subset of constituents pass in input
93
94	!! * Substitutions
95
96	!!----------------------------------------------------------------------
97	!! OPA 9.0 , LOCEAN-IPSL (2005)
98	!! or LIM 2.0 , UCL-LOCEAN-IPSL (2005)
99	!! or TOP 1.0 , LOCEAN-IPSL (2005)
100	!! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/module_example,v 1.3 2005/03/27 18:34:47 opalod Exp $
101	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
102	!!----------------------------------------------------------------------
103
104	CONTAINS
105
106	SUBROUTINE dia_harm_fast( kt )
107	!!----------------------------------------------------------------------
108	!! * ROUTINE harm_ana *
109	!!
110	!! ** Purpose : Harmonic analyser
111	!!
112	!! ** Method :
113	!!
114	!! ** Action : - first action (share memory array/varible modified
115	!! in this routine
116	!! - second action .....
117	!! - .....
118	!!
119	!! References :
120	!! Give references if exist otherwise suppress these lines
121	!!
122	!! History :
123	!! 9.0 ! 03-08 (Autor Names) Original code
124	!! ! 02-08 (Author names) brief description of modifications
125	!!----------------------------------------------------------------------
126	!! * Modules used
127
128	!! * arguments
129	INTEGER, INTENT( in ) :: &
130	kt ! describe it!!!
131
132	!! * local declarations
133	INTEGER :: ji, jk, jj ! dummy loop arguments
134	INTEGER :: jh, i1, i2, jgrid
135	INTEGER :: j2d, j3d
136	REAL(WP) :: sec2start,sec2start_old
137	CHARACTER (len=40) :: tmp_name
138	!!--------------------------------------------------------------------
139
140	!JT IF( nn_timing == 1 ) CALL timing_start( 'dia_harm_fast' )
141	IF( ln_timing ) CALL timing_start( 'dia_harm_fast' )
142	IF( kt == nit000 ) CALL harm_ana_init(kt) ! Initialization (first time-step only)
143
144
145	IF ( ln_diaharm_update_nodal_daily ) THEN
146	!IF (MOD(kt,nint(86400./rdt)) == 0) THEN
147
148
149	IF( nsec_day == NINT(0.5_wp * rdt) .OR. kt == nit000 ) THEN
150	DO jh = 1, nb_ana
151	!JT anau(jh) = 3.141579*utide ( ntide_sub(jh) )/180.
152	!JT anav(jh) = 3.141579*v0tide( ntide_sub(jh) )/180.
153	anau(jh) = utide ( ntide_sub(jh) )
154	anav(jh) = v0tide( ntide_sub(jh) )
155	anaf(jh) = ftide ( ntide_sub(jh) )
156	END DO
157
158	IF(lwp) WRITE(numout,*)
159	IF(lwp) WRITE(numout,*) 'harm_ana : update nodes?',ln_diaharm_update_nodal_daily
160	IF(lwp) WRITE(numout,*) 'harm_ana : date, time',ndastp, nhour, nminute
161
162	ENDIF
163	ENDIF
164
165
166	DO jh = 1, nb_ana
167
168	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_anau_utide'
169	IF( iom_use(TRIM(tmp_name)) ) THEN
170	! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(anau(jh) )
171	CALL iom_put( TRIM(tmp_name), anau(jh) )
172	!ELSE
173	! IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
174	ENDIF
175
176	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_anav_v0tide'
177	IF( iom_use(TRIM(tmp_name)) ) THEN
178	! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(anav(jh) )
179	CALL iom_put( TRIM(tmp_name), anav(jh) )
180	!ELSE
181	! IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
182	ENDIF
183
184	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_anaf_ftide'
185	IF( iom_use(TRIM(tmp_name)) ) THEN
186	! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(anaf(jh) )
187	CALL iom_put( TRIM(tmp_name), anaf(jh) )
188	!ELSE
189	! IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
190	ENDIF
191
192	END DO
193
194
195	IF ( ln_diaharm_fast .and. ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
196
197	! this bit done every time step
198	nhm=2*nb_ana+1
199	c(1) = 1.0
200
201	sec2start_old = nint( (fjulday-fjulday_startharm)*86400._wp )
202	sec2start = nsec_day - NINT(0.5_wp * rdt)
203	!IF(lwp) WRITE(numout,*) "ztime NEW", kt, sec2start, fjulday_startharm
204
205	DO jh=1,nb_ana
206	c(2jh ) = anaf(jh)cos( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
207	c(2jh+1) = anaf(jh)sin( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
208
209	IF(lwp) WRITE(numout,) 'diaharm_fast: analwave,',kt,tname(jh),sec2start,sec2start/3600.,sec2start_old,sec2start_old/3600,c(2jh),c(2*jh+1),om_tide(jh),anau(jh),anav(jh)
210	ENDDO
211
212	!IF(lwp) WRITE(numout,*) "c init", c, "c end", sec2start, om_tide(1), anau(1), anav(1),"end nodal"
213
214
215	! CUMULATE
216	DO ji=1,jpi ! loop lon
217	DO jj=1,jpj ! loop lat
218	DO jh=1,nhm ! loop harmonic
219
220	DO j2d=1,nvar_2d
221	IF ( m_posi_2d(j2d) .eq. 1 ) dd_cumul = c(jh) * sshn(ji,jj) * ssmask (ji,jj) ! analysis elevation
222	IF ( m_posi_2d(j2d) .eq. 2 ) dd_cumul = c(jh) * un_b(ji,jj) * ssumask(ji,jj) ! analysis depth average velocities
223	IF ( m_posi_2d(j2d) .eq. 3 ) dd_cumul = c(jh) * vn_b(ji,jj) * ssvmask(ji,jj)
224	!JT IF ( m_posi_2d(j2d) .eq. 4 ) dd_cumul = c(jh) * bfrua(ji,jj) * un(ji,jj,mbku(ji,jj)) * ssumask(ji,jj) ! analysis bottom friction
225	!JT IF ( m_posi_2d(j2d) .eq. 5 ) dd_cumul = c(jh) * bfrva(ji,jj) * vn(ji,jj,mbkv(ji,jj)) * ssvmask(ji,jj)
226	g_cumul_var2D(jh,ji,jj,j2d) = g_cumul_var2D(jh,ji,jj,j2d) + dd_cumul
227	ENDDO
228
229	DO j3d=1,nvar_3d
230	DO jk=1,jpkm1
231	IF ( m_posi_3d(j3d) .eq. 1 ) dd_cumul = c(jh) * rhd(ji,jj,jk) * tmask(ji,jj,jk)
232	IF ( m_posi_3d(j3d) .eq. 2 ) dd_cumul = c(jh) * ( un(ji,jj,jk)-un_b(ji,jj) ) * umask(ji,jj,jk)
233	IF ( m_posi_3d(j3d) .eq. 3 ) dd_cumul = c(jh) * ( vn(ji,jj,jk)-vn_b(ji,jj) ) * vmask(ji,jj,jk)
234	IF ( m_posi_3d(j3d) .eq. 4 ) dd_cumul = c(jh) * wn(ji,jj,jk) * wmask(ji,jj,jk)
235	g_cumul_var3D(jh,ji,jj,jk,j3d) = g_cumul_var3D(jh,ji,jj,jk,j3d) + dd_cumul
236	ENDDO
237	ENDDO
238
239	ENDDO ! end loop harmonic
240	ENDDO ! end loop lat
241	ENDDO ! end loop lon
242
243	! Compute nodal factor cumulative cross-product
244	DO i1=1,nhm
245	DO i2=1,nhm
246	cc(i1,i2)=cc(i1,i2)+c(i1)*c(i2)
247	ENDDO
248	ENDDO
249
250	! Output RESTART
251	IF( kt == nitrst ) THEN
252	CALL harm_rst_write(kt) ! Dump out data for a restarted run
253	ENDIF
254
255	! At End of run
256	IF ( kt == nitend ) THEN
257
258	IF(lwp) WRITE(numout,*)
259	IF(lwp) WRITE(numout,*) 'harm_ana : harmonic analysis of tides at end of run'
260	IF(lwp) WRITE(numout,*) '~~~~~~~~~'
261
262	IF( ln_diaharm_compute ) THEN
263
264	! INITIALISE TABLE TO 0
265	IF ( nvar_2d .gt. 0 ) THEN
266	g_cosamp2D = 0.0_wp
267	g_sinamp2D = 0.0_wp
268	ENDIF
269	IF ( nvar_3d .gt. 0 ) THEN
270	g_cosamp3D = 0.0_wp
271	g_sinamp3D = 0.0_wp
272	ENDIF
273
274	! FIRST OUTPUT 2D VARIABLES
275	DO jgrid=1,nvar_2d ! loop number of 2d variables (ssh, U2d, V2d, UVfric) to analyse harmonically
276	DO ji=1,jpi ! loop lon
277	DO jj=1,jpj ! loop lat
278	bt = 1.0_wp; bzz(:) = 0.0_wp
279	DO jh=1,nhm ! loop harmonic
280	bzz(jh) = g_cumul_var2D(jh,ji,jj,jgrid)
281	bt = bt*bzz(jh)
282	ENDDO
283	! Copy back original cumulated nodal factor
284	a(:,:) = cc(:,:)
285	! now do gaussian elimination of the system
286	! a * x = b
287	! the matrix x is (a0,a1,b1,a2,b2 ...)
288	! the matrix a and rhs b solved here for x
289	x=0.0_wp
290	IF(bt.ne.0.) THEN
291	CALL gelim( a, bzz, x, nhm )
292	! Backup output in variables
293	DO jh=1,nb_ana
294	g_cosamp2D(jh,ji,jj,jgrid) = x(jh*2 )
295	g_sinamp2D(jh,ji,jj,jgrid) = x(jh*2+1)
296	ENDDO
297	g_cosamp2D ( 0,ji,jj,jgrid) = x(1)
298	g_sinamp2D ( 0,ji,jj,jgrid) = 0.0_wp
299	ENDIF ! bt.ne.0.
300	ENDDO ! jj
301	ENDDO ! ji
302	ENDDO ! jgrid
303
304	! SECOND OUTPUT 3D VARIABLES
305	DO jgrid=1,nvar_3d ! loop number of 3d variables rho, U, V, W
306	DO jk=1,jpkm1 ! loop over vertical level
307	DO ji=1,jpi ! loop over lon
308	DO jj=1,jpj ! loop over lat
309	bt = 1.0_wp; bzz(:) = 0.0_wp
310	DO jh=1,nhm
311	bzz(jh) = g_cumul_var3D(jh,ji,jj,jk,jgrid)
312	bt = bt*bzz(jh)
313	ENDDO
314	! Copy back original cumulated nodal factor
315	a(:,:) = cc(:,:)
316	! now do gaussian elimination of the system
317	! a * x = b
318	! the matrix x is (a0,a1,b1,a2,b2 ...)
319	! the matrix a and rhs b solved here for x
320	x=0.0_wp
321	IF(bt.ne.0.) THEN
322	CALL gelim( a, bzz, x, nhm )
323	! Backup output in variables
324	DO jh=1,nb_ana
325	g_cosamp3D(jh,ji,jj,jk,jgrid) = x(jh*2 )
326	g_sinamp3D(jh,ji,jj,jk,jgrid) = x(jh*2+1)
327	ENDDO
328	g_cosamp3D ( 0,ji,jj,jk,jgrid) = x(1)
329	g_sinamp3D ( 0,ji,jj,jk,jgrid) = 0.0_wp
330	ENDIF ! bt.ne.0.
331	ENDDO ! jj
332	ENDDO ! ji
333	ENDDO ! jk
334	ENDDO ! jgrid
335
336	CALL harm_ana_out ! output analysis (last time step)
337
338	ELSE ! ln_harmana_compute = False
339	IF(lwp) WRITE(numout,*) " Skipping Computing harmonics at last step"
340
341	ENDIF ! ln_harmana_compute
342	ENDIF ! kt == nitend
343
344	ENDIF
345
346	!JT IF( nn_timing == 1 ) CALL timing_stop( 'dia_harm_fast' )
347	IF( ln_timing ) CALL timing_stop( 'dia_harm_fast' )
348
349	END SUBROUTINE dia_harm_fast
350
351	SUBROUTINE harm_ana_init( kt ) !JT
352	!!----------------------------------------------------------------------
353	!! * ROUTINE harm_ana_init *
354	!!
355	!! ** Purpose : initialization of ....
356	!!
357	!! ** Method : blah blah blah ...
358	!!
359	!! ** input : Namlist namexa
360	!!
361	!! ** Action : ...
362	!!
363	!! history :
364	!! 9.0 ! 03-08 (Autor Names) Original code
365	!!----------------------------------------------------------------------
366	!! * local declarations
367
368	INTEGER, INTENT(in) :: kt ! ocean time-step !JT
369	!!
370	INTEGER :: ji, jk, jh ! dummy loop indices
371	INTEGER :: ios ! Local integer output status for namelist read
372	INTEGER :: k2d, k3d ! dummy number of analysis
373
374
375
376
377	NAMELIST/nam_diaharm_fast/ ln_diaharm_fast, ln_diaharm_store, ln_diaharm_compute, ln_diaharm_read_restart, ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d, &
378	& tname,ln_diaharm_multiyear,nn_diaharm_multiyear,ln_diaharm_update_nodal_daily
379	!!----------------------------------------------------------------------
380	!JT
381	ln_diaharm_fast = .FALSE.
382	ln_diaharm_multiyear = .FALSE.
383	nn_diaharm_multiyear = 20
384	!JT
385	lk_diaharm_2D = .TRUE. ! to run 2d
386	lk_diaharm_3D = .TRUE. ! to run 3d
387
388	ln_diaharm_store = .TRUE.
389
390	IF(lwp) WRITE(numout,*)
391	IF(lwp) WRITE(numout,*) 'harm_init : initialization of harmonic analysis of tides'
392	IF(lwp) WRITE(numout,*) '~~~~~~~~~'
393
394	! GET NAMELIST DETAILS
395	REWIND( numnam_ref ) ! Namelist nam_diaharm_fast in reference namelist : Tidal harmonic analysis
396	READ ( numnam_ref, nam_diaharm_fast, IOSTAT = ios, ERR = 901)
397	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in reference namelist' )
398
399	REWIND( numnam_cfg ) ! Namelist nam_diaharm_fast in configuration namelist : Tidal harmonic analysis
400	READ ( numnam_cfg, nam_diaharm_fast, IOSTAT = ios, ERR = 902 )
401	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in configuration namelist' )
402	IF(lwm) WRITE ( numond, nam_diaharm_fast )
403
404
405	!
406	IF(lwp) THEN
407	WRITE(numout,*) 'Tidal diagnostics_fast'
408	WRITE(numout,*) ' Fast Harmonic analysis ?: ln_diaharm_fast= ', ln_diaharm_fast
409	WRITE(numout,*) ' Store output in restart?: ln_diaharm_store= ', ln_diaharm_store
410	WRITE(numout,*) ' Compute analysis?: ln_diaharm_compute= ', ln_diaharm_compute
411	WRITE(numout,*) ' Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
412	WRITE(numout,*) ' SSH harmonic analysis: ln_ana_ssh = ', ln_ana_ssh
413	WRITE(numout,*) ' Barotropic Velocities harmonic analysis: ln_ana_uvbar = ', ln_ana_uvbar
414	WRITE(numout,*) ' Bed Friction for harmonic analysis (not implemented): ln_ana_bfric = ', ln_ana_bfric
415	WRITE(numout,*) ' Density harmonic analysis: ln_ana_rho = ', ln_ana_rho
416	WRITE(numout,*) ' 3D velocities harmonic analysis: ln_ana_uv3d = ', ln_ana_uv3d
417	WRITE(numout,*) ' Vertical Velocities harmonic analysis: ln_ana_w3d = ', ln_ana_w3d
418	WRITE(numout,*) ' Names of harmonics: tname = ', tname
419	WRITE(numout,*) ' Max number of harmonics: jpmax_harmo = ', jpmax_harmo
420	WRITE(numout,*) ' Number of Harmonics: nb_harmo = ', nb_harmo
421	WRITE(numout,*) ' Multi-year harmonic analysis: ln_diaharm_multiyear = ', ln_diaharm_multiyear
422	WRITE(numout,*) ' Multi-year harmonic analysis - number of years: nn_diaharm_multiyear = ', nn_diaharm_multiyear
423	WRITE(numout,*) ' Multi-year harmonic analysis - number of years: ln_diaharm_update_nodal_daily = ', ln_diaharm_update_nodal_daily
424	WRITE(numout,*) ' Number of Harmonics: nyear, nmonth = ', nyear, nmonth
425
426	ENDIF
427	! JT
428
429
430	IF ( ln_diaharm_multiyear ) THEN
431	ln_diaharm_store = .True.
432	IF(lwp) WRITE(numout,*) ' Multi-year harmonic analysis ', nyear,nn_diaharm_multiyear,nmonth
433	IF ((mod(nyear,nn_diaharm_multiyear) == 0) .AND. ( nmonth == 1 )) THEN ! Jan, year = 1980,2000,2020,2040, restart tidal calculation
434	ln_diaharm_read_restart = .FALSE.
435	IF(lwp) WRITE(numout,*) ' Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
436	ELSE
437	ln_diaharm_read_restart = .TRUE.
438	IF(lwp) WRITE(numout,*) ' Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
439	ENDIF
440
441
442
443	IF(lwp) WRITE(numout,*) ' Multi-year harmonic analysis ', nyear,nn_diaharm_multiyear,nmonth
444	IF ((mod(nyear,nn_diaharm_multiyear) == (nn_diaharm_multiyear - 1)) .AND. ( nmonth == 12 )) THEN ! Dec year = 1999,2019,2039,2040, restart tidal calculation
445	ln_diaharm_compute = .TRUE.
446	IF(lwp) WRITE(numout,*) ' Compute analysis?: ln_diaharm_compute= ', ln_diaharm_compute
447	ELSE
448	ln_diaharm_compute = .FALSE.
449	IF(lwp) WRITE(numout,*) ' Compute analysis?: ln_diaharm_compute= ', ln_diaharm_compute
450	ENDIF
451
452	ENDIF
453	IF ( kt < 10 ) THEN
454	ln_diaharm_read_restart = .FALSE.
455	IF(lwp) WRITE(numout,*) ' kt = ',kt
456	IF(lwp) WRITE(numout,*) ' kt < 10, so setting ln_diaharm_read_restart to .FALSE.'
457	IF(lwp) WRITE(numout,*) ' Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
458	ENDIF
459
460	! JT
461
462
463	! GET NUMBER OF HARMONIC TO ANALYSE - from diaharm.F90
464	nb_ana = 0
465	DO jk=1,jpmax_harmo
466	DO ji=1,nb_harmo
467	IF(TRIM(tname(jk)) == TRIM(Wave( ntide(ji) )%cname_tide) ) THEN
468	nb_ana=nb_ana+1
469	ENDIF
470	END DO
471	END DO
472	!
473	IF(lwp) THEN
474	WRITE(numout,*) ' Namelist nam_diaharm_fast'
475	WRITE(numout,*) ' nb_ana = ', nb_ana
476	CALL flush(numout)
477	ENDIF
478	!
479	IF (nb_ana > nharm_max) THEN
480	IF(lwp) WRITE(numout,*) ' E R R O R harm_ana : nb_ana must be lower than nharm_max, stop'
481	IF(lwp) WRITE(numout,*) ' nharm_max = ', nharm_max
482	nstop = nstop + 1
483	ENDIF
484
485
486	ALLOCATE(ntide_all(nb_ana))
487	ALLOCATE(ntide_sub(nb_ana))
488
489	DO jk=1,nb_ana
490	DO ji=1,nb_harmo
491	IF (TRIM(tname(jk)) .eq. Wave( ntide(ji) )%cname_tide ) THEN
492	ntide_sub(jk) = ji
493	ntide_all(jk) = ntide(ji)
494	EXIT
495	END IF
496	END DO
497	END DO
498
499	ALLOCATE( anau(nb_ana) )
500	ALLOCATE( anav(nb_ana) )
501	ALLOCATE( anaf(nb_ana) )
502
503
504	IF( ln_diaharm_fast ) THEN
505
506	! SEARCH HOW MANY VARIABLES 2D AND 3D TO PROCESS
507	nvar_2d = 0; nvar_3d = 0
508	IF ( ln_ana_ssh ) nvar_2d = nvar_2d + 1 ! analysis elevation
509	IF ( ln_ana_uvbar ) nvar_2d = nvar_2d + 2 ! analysis depth-averaged velocity
510	IF ( ln_ana_bfric ) nvar_2d = nvar_2d + 2 ! analysis bottom friction
511
512	IF ( ln_ana_rho ) nvar_3d = nvar_3d + 1 ! analysis density
513	IF ( ln_ana_uv3d ) nvar_3d = nvar_3d + 2 ! analysis 3D horizontal velocities
514	IF ( ln_ana_w3d ) nvar_3d = nvar_3d + 1 ! analysis 3D vertical velocity
515
516	! CHECK IF SOMETHING TO RUN
517	IF ( nvar_2d .eq. 0 ) lk_diaharm_2D = .FALSE. ! no 2d to run
518	IF ( nvar_3d .eq. 0 ) lk_diaharm_3D = .FALSE. ! no 3d to run
519	! IF ( nvar_2d .gt. 0 .and. nvar_3d .gt. 0 ) lk_diaharm_fast = .FALSE.
520	! IF ( .NOT. ln_diaharm_store ) lk_diaharm_fast = .FALSE.
521
522	IF ( ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
523
524	! DO ALLOCATIONS
525	IF ( lk_diaharm_2D ) THEN
526	ALLOCATE( g_cumul_var2D(nb_ana*2+1,jpi,jpj, nvar_2d) )
527	ALLOCATE( g_cosamp2D( 0:nb_ana*2+1,jpi,jpj, nvar_2d) )
528	ALLOCATE( g_sinamp2D( 0:nb_ana*2+1,jpi,jpj, nvar_2d) )
529	ALLOCATE( g_out2D (jpi,jpj) )
530	ALLOCATE( h_out2D (jpi,jpj) )
531	ALLOCATE( m_posi_2d( nvar_2d ) ); m_posi_2d(:)=0
532	ENDIF
533
534	IF ( lk_diaharm_3D ) THEN
535	ALLOCATE( g_cumul_var3D(nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
536	ALLOCATE( g_cosamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
537	ALLOCATE( g_sinamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
538	ALLOCATE( g_out3D (jpi,jpj,jpk) )
539	ALLOCATE( h_out3D (jpi,jpj,jpk) )
540	ALLOCATE( m_posi_3d( nvar_3d ) ); m_posi_3d(:)=0
541	ENDIF
542
543	ALLOCATE( cc(nb_ana2+1,nb_ana2+1) )
544	ALLOCATE( a (nb_ana2+1,nb_ana2+1) )
545	ALLOCATE( bzz(nb_ana*2+1) )
546	ALLOCATE( x (nb_ana*2+1) )
547	ALLOCATE( c (nb_ana*2+1) )
548	! END ALLOCATE
549
550	! STORE INDEX OF WHAT TO PRODUCE DEPENDING ON ACTIVATED LOGICAL
551	! MAKES THINGS EASIER AND FASTER LATER
552	! !!! UGLY !!!
553	jh = 1; k2d = 0;
554	IF ( ln_ana_ssh ) THEN
555	k2d = k2d + 1; m_posi_2d(k2d) = jh
556	IF(lwp) WRITE(numout,*) " - ssh harmonic analysis activated (ln_ana_ssh)"
557	ENDIF
558	jh = jh + 1
559	IF ( ln_ana_uvbar ) THEN
560	k2d = k2d + 1; m_posi_2d(k2d) = jh
561	jh = jh + 1
562	k2d = k2d + 1; m_posi_2d(k2d) = jh
563	IF(lwp) WRITE(numout,*) " - barotropic currents harmonic analysis activated (ln_ana_uvbar)"
564	ELSE
565	jh = jh + 1
566	ENDIF
567	jh = jh + 1
568	IF ( ln_ana_bfric ) THEN
569	k2d = k2d + 1; m_posi_2d(k2d) = jh
570	jh = jh + 1;
571	k2d = k2d + 1; m_posi_2d(k2d) = jh
572	IF(lwp) WRITE(numout,*) " - bottom friction harmonic analysis activated (ln_ana_vbfr)"
573	ELSE
574	jh = jh + 1
575	ENDIF
576
577	! and for 3D
578	jh = 1; k3d = 0;
579	IF ( ln_ana_rho ) THEN
580	k3d = k3d + 1; m_posi_3d(k3d) = jh
581	IF(lwp) WRITE(numout,*) " - 3D density harmonic analysis activated (ln_ana_rho)"
582	ENDIF
583	jh = jh + 1
584	IF ( ln_ana_uv3d ) THEN
585	k3d = k3d + 1; m_posi_3d(k3d) = jh
586	jh = jh + 1
587	k3d = k3d + 1; m_posi_3d(k3d) = jh
588	IF(lwp) WRITE(numout,*) " - 3D horizontal currents harmonic analysis activated (ln_ana_uv3d)"
589	ELSE
590	jh = jh + 1
591	ENDIF
592	jh = jh + 1
593	IF ( ln_ana_w3d ) THEN
594	k3d = k3d + 1; m_posi_3d(k3d) = jh
595	IF(lwp) WRITE(numout,*) " - 3D vertical currents harmonic analysis activated (ln_ana_w3d)"
596	ENDIF
597
598	! SELECT AND STORE FREQUENCIES
599	IF(lwp) WRITE(numout,*) 'Analysed frequency : ',nb_ana ,'Frequency '
600	DO jh=1,nb_ana
601	om_tide(jh) = omega_tide( ntide_sub(jh) )
602	IF(lwp) WRITE(numout,) ' - ',tname(jh),' ',om_tide(jh), (2rpi/3600.)/om_tide(jh),"hr"
603	ENDDO
604
605	! READ RESTART IF
606	IF ( ln_diaharm_read_restart ) THEN
607	IF (lwp) WRITE(numout,*) "Reading previous harmonic data from previous run. kt = ",kt
608	! Need to read in bssh bz, cc anau anav and anaf
609	call harm_rst_read ! This reads in from the previous day
610	! Currrently the data in in assci format
611	ELSE
612
613	IF (lwp) WRITE(numout,*) "Starting harmonic analysis from Fresh. kt = ",kt
614
615	IF ( lk_diaharm_2D ) g_cumul_var2D(:,:,:,: ) = 0.0_wp
616	IF ( lk_diaharm_3D ) g_cumul_var3D(:,:,:,:,:) = 0.0_wp
617	cc = 0.0_wp
618	a (:,:) = 0.0_wp ! NB
619	bzz (:) = 0.0_wp
620	x (:) = 0.0_wp
621	c (:) = 0.0_wp
622	anau (:) = 0.0_wp
623	anav (:) = 0.0_wp
624	anaf (:) = 0.0_wp
625
626	DO jh = 1, nb_ana
627	anau(jh) = utide ( ntide_sub(jh) )
628	anav(jh) = v0tide( ntide_sub(jh) )
629	anaf(jh) = ftide ( ntide_sub(jh) )
630	END DO
631
632	fjulday_startharm=fjulday !Set this at very start and store
633	!JT this is a mistake - only works on daily cycles, should use fjulnsec_dayday
634
635	IF (lwp) THEN
636	WRITE(numout,*) '--------------------------'
637	WRITE(numout,*) ' - Output anaf for check'
638	WRITE(numout,*) 'ANA F', anaf
639	WRITE(numout,*) 'ANA U', anau
640	WRITE(numout,*) 'ANA V', anav
641	WRITE(numout,*) 'fjulday',fjulday
642	WRITE(numout,*) 'fjulday_startharm',fjulday_startharm
643	WRITE(numout,*) 'nsec_day',nsec_day
644	WRITE(numout,*) 'kt',kt
645	WRITE(numout,*) '--------------------------'
646	ENDIF
647
648	ENDIF
649
650	ELSE
651
652	IF (lwp) WRITE(numout,*) "No variable setup for harmonic analysis"
653
654	ENDIF
655	ENDIF
656
657	END SUBROUTINE harm_ana_init
658	!
659	SUBROUTINE gelim (a,b,x,n)
660	!!----------------------------------------------------------------------
661	!! * ROUTINE harm_ana *
662	!!
663	!! ** Purpose : Guassian elimination
664	!!
665	!!
666	!! ** Action : - first action (share memory array/varible modified
667	!! in this routine
668	!! - second action .....
669	!! - .....
670	!!
671	!! References :
672	!! Give references if exist otherwise suppress these lines
673	!!
674	!! History :
675	implicit none
676	!
677	integer :: n
678	REAL(WP) :: b(nb_ana2+1), a(nb_ana2+1,nb_ana*2+1)
679	REAL(WP) :: x(nb_ana*2+1)
680	INTEGER :: row,col,prow,pivrow,rrow
681	REAL(WP) :: atemp
682	REAL(WP) :: pivot
683	REAL(WP) :: m
684
685	do row=1,n-1
686	pivrow=row
687	pivot=a(row,n-row+1)
688	do prow=row+1,n
689	if (abs(a(prow,n-row+1)).gt.abs(pivot) ) then
690	pivot=a(prow,n-row+1)
691	pivrow=prow
692	endif
693	enddo
694	! swap row and prow
695	if ( pivrow .ne. row ) then
696	atemp=b(pivrow)
697	b(pivrow)=b(row)
698	b(row)=atemp
699	do col=1,n
700	atemp=a(pivrow,col)
701	a(pivrow,col)=a(row,col)
702	a(row,col)=atemp
703	enddo
704	endif
705
706	do rrow=row+1,n
707	if (a(row,row).ne.0) then
708
709	m=-a(rrow,n-row+1)/a(row,n-row+1)
710	do col=1,n
711	a(rrow,col)=m*a(row,col)+a(rrow,col)
712	enddo
713	b(rrow)=m*b(row)+b(rrow)
714	endif
715	enddo
716	enddo
717	! back substitution now
718
719	x(1)=b(n)/a(n,1)
720	do row=n-1,1,-1
721	x(n-row+1)=b(row)
722	do col=1,(n-row)
723	x(n-row+1)=(x(n-row+1)-a(row,col)*x(col))
724	enddo
725
726	x(n-row+1)=(x(n-row+1)/a(row,(n-row)+1))
727	enddo
728
729	return
730	END SUBROUTINE gelim
731
732	SUBROUTINE harm_ana_out
733	!!----------------------------------------------------------------------
734	!! * ROUTINE harm_ana_init *
735	!!
736	!! ** Purpose : initialization of ....
737	!!
738	!! ** Method : blah blah blah ...
739	!!
740	!! ** input : Namlist namexa
741	!!
742	!! ** Action : ...
743	!!
744	!! history :
745	!! 9.0 ! 03-08 (Autor Names) Original code
746	!!----------------------------------------------------------------------
747	USE dianam ! build name of file (routine)
748
749	!! * local declarations
750	INTEGER :: ji, jj, jk, jgrid, jh ! dummy loop indices
751	! INTEGER :: nh_T
752	! INTEGER :: nid_harm
753	! CHARACTER (len=40) :: clhstnamt, clop1, clop2 ! temporary names
754	! CHARACTER (len=40) :: clhstnamu, clhstnamv ! temporary names
755	CHARACTER (len=40) :: suffix
756	CHARACTER (len=40) :: tmp_name
757	! REAL(wp) :: zsto1, zsto2, zout, zmax, zjulian, zdt, zmdi ! temporary scalars
758
759	do jgrid=1,nvar_2d
760	do jh=1,nb_ana
761	h_out2D = 0.0
762	g_out2D = 0.0
763	do jj=1,nlcj
764	do ji=1,nlci
765	cca=g_cosamp2D(jh,ji,jj,jgrid)
766	ssa=g_sinamp2D(jh,ji,jj,jgrid)
767	h_out2D(ji,jj)=sqrt(cca2+ssa2)
768	IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
769	g_out2D(ji,jj)= 0.0_wp
770	ELSE
771	g_out2D(ji,jj)=(180.0/rpi)*atan2(ssa,cca)
772	ENDIF
773	IF (h_out2D(ji,jj).ne.0) THEN
774	h_out2D(ji,jj)=h_out2D(ji,jj)/anaf(jh)
775	ENDIF
776	IF (g_out2D(ji,jj).ne.0) THEN !Correct and take modulus
777	g_out2D(ji,jj) = g_out2D(ji,jj) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
778	if (g_out2D(ji,jj).gt.360.0) then
779	g_out2D(ji,jj)=g_out2D(ji,jj)-360.0
780	else if (g_out2D(ji,jj).lt.0.0) then
781	g_out2D(ji,jj)=g_out2D(ji,jj)+360.0
782	endif
783	ENDIF
784	enddo
785	enddo
786	!
787	! NETCDF OUTPUT
788	suffix = TRIM( m_varName2d( m_posi_2d(jgrid) ) )
789	IF(lwp) WRITE(numout,*) "harm_ana_out", suffix
790
791	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix)
792	IF( iom_use(TRIM(tmp_name)) ) THEN
793	IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out2D)
794	IF(lwp) WRITE(numout,) "harm_ana_out_names", tmp_name,tname(jh),' ',om_tide(jh), (2rpi/3600.)/om_tide(jh),"hr"
795	CALL iom_put( TRIM(tmp_name), h_out2D(:,:) )
796	ELSE
797	IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
798	ENDIF
799
800	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix)
801	IF( iom_use(TRIM(tmp_name)) ) THEN
802	IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out2D)
803	CALL iom_put( TRIM(tmp_name), g_out2D(:,:) )
804	ELSE
805	IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
806	ENDIF
807
808	enddo
809
810	!tmp_name=TRIM(suffix)//'_tidal_anal_offset'
811	!IF( iom_use(TRIM(tmp_name)) ) THEN
812	! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name)
813	! CALL iom_put( TRIM(tmp_name), g_cosamp2D( 0,:,:,jgrid))
814	!ELSE
815	! IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
816	!ENDIF
817
818	enddo
819	!
820	! DO THE SAME FOR 3D VARIABLES
821	!
822	do jgrid=1,nvar_3d
823	do jh=1,nb_ana
824	h_out3D = 0.0
825	g_out3D = 0.0
826	DO jk=1,jpkm1
827	do jj=1,nlcj
828	do ji=1,nlci
829	cca=g_cosamp3D(jh,ji,jj,jk,jgrid)
830	ssa=g_sinamp3D(jh,ji,jj,jk,jgrid)
831	h_out3D(ji,jj,jk)=sqrt(cca2+ssa2)
832	IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
833	g_out3D(ji,jj,jk) = 0.0_wp
834	ELSE
835	g_out3D(ji,jj,jk) = (180.0/rpi)*atan2(ssa,cca)
836	ENDIF
837	IF (h_out3D(ji,jj,jk).ne.0) THEN
838	h_out3D(ji,jj,jk) = h_out3D(ji,jj,jk)/anaf(jh)
839	ENDIF
840	IF (g_out3D(ji,jj,jk).ne.0) THEN !Correct and take modulus
841	g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
842	if (g_out3D(ji,jj,jk).gt.360.0) then
843	g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)-360.0
844	else if (g_out3D(ji,jj,jk).lt.0.0) then
845	g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)+360.0
846	endif
847	ENDIF
848	enddo ! ji
849	enddo ! jj
850	ENDDO ! jk
851	!
852	! NETCDF OUTPUT
853	suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) )
854	IF(lwp) WRITE(numout,*) "harm_ana_out", suffix
855
856	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix)
857	IF( iom_use(TRIM(tmp_name)) ) THEN
858	IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out3D)
859	CALL iom_put( TRIM(tmp_name), h_out3D(:,:,:) )
860	ELSE
861	IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
862	ENDIF
863
864	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix)
865	IF( iom_use(TRIM(tmp_name)) ) THEN
866	IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out3D)
867	CALL iom_put(tmp_name, g_out3D(:,:,:) )
868	ELSE
869	IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
870	ENDIF
871
872	enddo ! jh
873
874	!tmp_name=TRIM(suffix)//'_tidal_anal_offset'
875	!IF( iom_use(TRIM(tmp_name)) ) THEN
876	! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name)
877	! CALL iom_put( TRIM(tmp_name), g_cosamp3D( 0,:,:,:,jgrid))
878	!ELSE
879	! IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
880	!ENDIF
881
882	enddo ! jgrid
883	!
884	END SUBROUTINE harm_ana_out
885	!
886	SUBROUTINE harm_rst_write(kt)
887	!!----------------------------------------------------------------------
888	!! * ROUTINE harm_ana_init *
889	!!
890	!! ** Purpose : To write out cummulated Tidal Harmomnic data to file for
891	!! restarting
892	!!
893	!! ** Method : restart files will be dated by default
894	!!
895	!! ** input :
896	!!
897	!! ** Action : ...
898	!!
899	!! history :
900	!! 0.0 ! 01-16 (Enda O'Dea) Original code
901	!! ASSUMES dated file for rose , can change later to be more generic
902	!!----------------------------------------------------------------------
903	INTEGER, INTENT(in) :: kt ! ocean time-step
904	!!
905	INTEGER :: jh, j2d, j3d
906	CHARACTER(LEN=20) :: clkt ! ocean time-step define as a character
907	CHARACTER(LEN=50) :: clname ! ocean output restart file name
908	CHARACTER(LEN=150) :: clpath ! full path to ocean output restart file
909	CHARACTER(LEN=250) :: clfinal ! full name
910
911	!restart file
912	DO j2d=1,nvar_2d
913	CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
914	DO jh=1,nb_ana
915	CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2 , :, :, j2d ) )
916	CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) )
917	ENDDO
918	ENDDO
919
920	DO j3d=1,nvar_3d
921	!JT CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
922	CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName3d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
923	DO jh=1,nb_ana
924	CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2 , :, :, :, j3d ) )
925	CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) )
926	ENDDO
927	ENDDO
928
929	IF(lwp) THEN
930	IF( kt > 999999999 ) THEN ; WRITE(clkt, * ) kt
931	ELSE ; WRITE(clkt, '(i8.8)') kt
932	ENDIF
933	clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
934	clpath = TRIM(cn_ocerst_outdir)
935	IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
936	IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for writing: ',TRIM(clpath)//clname
937
938	WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
939	OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
940	WRITE(66) cc
941	WRITE(66) anau
942	WRITE(66) anav
943	WRITE(66) anaf
944	WRITE(66) fjulday_startharm
945	CLOSE(66)
946	WRITE(numout,*) '----------------------------'
947	WRITE(numout,*) ' harm_rst_write: DONE '
948	WRITE(numout,*) cc
949	WRITE(numout,*) anaf
950	WRITE(numout,*) anau
951	WRITE(numout,*) anav
952	WRITE(numout,*) fjulday_startharm
953	WRITE(numout,*) '----------------------------'
954	ENDIF
955
956	END SUBROUTINE harm_rst_write
957
958	SUBROUTINE harm_rst_read
959	!!----------------------------------------------------------------------
960	!! * ROUTINE harm_ana_init *
961	!!
962	!! ** Purpose : To read in cummulated Tidal Harmomnic data to file for
963	!! restarting
964	!!
965	!! ** Method :
966	!!
967	!! ** input :
968	!!
969	!! ** Action : ...
970	!!
971	!! history :
972	!! 0.0 ! 01-16 (Enda O'Dea) Original code
973	!! ASSUMES dated file for rose , can change later to be more generic
974	!!----------------------------------------------------------------------
975	CHARACTER(LEN=20) :: clkt ! ocean time-step define as a character
976	CHARACTER(LEN=50) :: clname ! ocean output restart file name
977	CHARACTER(LEN=150) :: clpath ! full path to ocean output restart file
978	CHARACTER(LEN=250) :: clfinal ! full name
979	INTEGER :: jh, j2d, j3d
980
981	IF( nit000 > 999999999 ) THEN ; WRITE(clkt, * ) nit000-1
982	ELSE ; WRITE(clkt, '(i8.8)') nit000-1
983	ENDIF
984	clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
985	clpath = TRIM(cn_ocerst_outdir)
986	IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
987
988	IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for reading: ',TRIM(clpath)//clname
989
990	DO j2d=1,nvar_2d
991	CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
992	IF(lwp) WRITE(numout,*) "2D", j2d, m_posi_2d(j2d), m_varName2d( m_posi_2d(j2d) )
993	DO jh=1,nb_ana
994	CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2 , :, :, j2d ) )
995	CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) )
996	ENDDO
997	ENDDO
998
999	DO j3d=1,nvar_3d
1000	!JT CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
1001	CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName3d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
1002	IF(lwp) WRITE(numout,*) "3D", j3d, m_posi_3d(j3d), m_varName3d( m_posi_3d(j3d) )
1003
1004	DO jh=1,nb_ana
1005	CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2 , :, :, :, j3d ) )
1006	CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) )
1007	ENDDO
1008	ENDDO
1009
1010	WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
1011	OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
1012	READ(66) cc
1013	READ(66) anau
1014	READ(66) anav
1015	READ(66) anaf
1016	READ(66) fjulday_startharm
1017	CLOSE(66)
1018
1019	IF(lwp) THEN
1020	WRITE(numout,*) '----------------------------'
1021	WRITE(numout,*) ' Checking anaf is correct'
1022	WRITE(numout,*) cc
1023	WRITE(numout,*) anaf
1024	WRITE(numout,*) fjulday_startharm
1025	WRITE(numout,*) '----------------------------'
1026	ENDIF
1027
1028	END SUBROUTINE harm_rst_read
1029
1030	!!======================================================================
1031	#else
1032	!!---------------------------------------------------------------------------------
1033	!! Dummy module NO harmonic Analysis
1034	!!---------------------------------------------------------------------------------
1035	LOGICAL, PUBLIC, PARAMETER :: lk_diaharm_fast = .FALSE. ! to be run or not
1036
1037	CONTAINS
1038	SUBROUTINE harm_rst_write(kt) ! Dummy routine
1039	END SUBROUTINE harm_rst_write
1040	SUBROUTINE harm_rst_read ! Dummy routine
1041	END SUBROUTINE harm_rst_read
1042	SUBROUTINE harm_ana_out ! Dummy routine
1043	END SUBROUTINE harm_ana_out
1044	SUBROUTINE harm_ana_init
1045	END SUBROUTINE harm_ana_init
1046	SUBROUTINE harm_ana( kt )
1047	!--- NB : end call not properly written
1048	END SUBROUTINE harm_ana
1049	! END SUBROUTINE harm_ana_init
1050	!--- END NB
1051	SUBROUTINE gelim (a,b,x,n)
1052	!--- NB : end call not properly written
1053	END SUBROUTINE gelim
1054	! END SUBROUTINE gelim (a,b,x,n)
1055	!--- END NB
1056	#endif
1057
1058	END MODULE diaharm_fast

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA/diaharm_fast.F90 @ 15500

Download in other formats: