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.

diaharm_fast.F90 in NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA – NEMO

Context Navigation

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

Last change on this file since 15547 was 15547, checked in by hadjt, 3 years ago
Moving ftide, utide, v0tide output to sbctides, adding tide_t output to diaharm_fast
File size: 44.9 KB

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

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

Download in other formats: