Context Navigation

diaharm_fast.F90 @ 15491

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

DIA/diaharm_fast.F90

Harmonic analysis with restarts.

Working, but perhaps needs additional work

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

Download in other formats: