Context Navigation

diaharm_fast.F90 @ 15608

Last change on this file since 15608 was 15608, checked in by hadjt, 2 years ago

DIA/diaharm_fast.F90 working

bug fixed diaharm_fast.F90

amp_u3d(jh,ji,jj,jk) = h_out3D(ji,jj,jk) was outside the harmonics loop, so only ran for M2
this also fixed the vertical stripes in the 3d output fields

Code also tidied up, and a namelist verbosity switch was added

File size: 87.4 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	LOGICAL, PUBLIC :: ln_diaharm_postproc_vel
85	LOGICAL, PUBLIC :: ln_diaharm_verbose
86
87
88	!JT
89	LOGICAL, PUBLIC :: ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d
90	INTEGER :: nb_ana ! Number of harmonics to analyse
91	CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: tname ! Names of tidal constituents ('M2', 'K1',...)
92	INTEGER , ALLOCATABLE, DIMENSION(:) :: ntide_all ! INDEX within the full set of constituents (tide.h90)
93	INTEGER , ALLOCATABLE, DIMENSION(:) :: ntide_sub ! INDEX within the subset of constituents pass in input
94
95	!! * Substitutions
96
97	!!----------------------------------------------------------------------
98	!! OPA 9.0 , LOCEAN-IPSL (2005)
99	!! or LIM 2.0 , UCL-LOCEAN-IPSL (2005)
100	!! or TOP 1.0 , LOCEAN-IPSL (2005)
101	!! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/module_example,v 1.3 2005/03/27 18:34:47 opalod Exp $
102	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
103	!!----------------------------------------------------------------------
104
105	CONTAINS
106
107	SUBROUTINE dia_harm_fast( kt )
108	!!----------------------------------------------------------------------
109	!! * ROUTINE harm_ana *
110	!!
111	!! ** Purpose : Harmonic analyser
112	!!
113	!! ** Method :
114	!!
115	!! ** Action : - first action (share memory array/varible modified
116	!! in this routine
117	!! - second action .....
118	!! - .....
119	!!
120	!! References :
121	!! Give references if exist otherwise suppress these lines
122	!!
123	!! History :
124	!! 9.0 ! 03-08 (Autor Names) Original code
125	!! ! 02-08 (Author names) brief description of modifications
126	!!----------------------------------------------------------------------
127	!! * Modules used
128
129	!! * arguments
130	INTEGER, INTENT( in ) :: &
131	kt ! describe it!!!
132
133	!! * local declarations
134	INTEGER :: ji, jk, jj ! dummy loop arguments
135	INTEGER :: jh, i1, i2, jgrid
136	INTEGER :: j2d, j3d
137	REAL(WP) :: sec2start,sec2start_old
138	CHARACTER (len=40) :: tmp_name
139	!!--------------------------------------------------------------------
140
141	!JT IF( nn_timing == 1 ) CALL timing_start( 'dia_harm_fast' )
142	IF( ln_timing ) CALL timing_start( 'dia_harm_fast' )
143	IF( kt == nit000 ) CALL harm_ana_init(kt) ! Initialization (first time-step only)
144
145
146	IF ( ln_diaharm_update_nodal_daily ) THEN
147	!IF (MOD(kt,nint(86400./rdt)) == 0) THEN
148
149
150	IF( nsec_day == NINT(0.5_wp * rdt) .OR. kt == nit000 ) THEN
151	DO jh = 1, nb_ana
152	!JT anau(jh) = 3.141579*utide ( ntide_sub(jh) )/180.
153	!JT anav(jh) = 3.141579*v0tide( ntide_sub(jh) )/180.
154	anau(jh) = utide ( ntide_sub(jh) )
155	anav(jh) = v0tide( ntide_sub(jh) )
156	anaf(jh) = ftide ( ntide_sub(jh) )
157	END DO
158
159	IF(lwp) WRITE(numout,*)
160	IF(lwp) WRITE(numout,*) 'harm_ana : update nodes?',ln_diaharm_update_nodal_daily
161	IF(lwp) WRITE(numout,*) 'harm_ana : date, time',ndastp, nhour, nminute
162
163	ENDIF
164	ENDIF
165
166	!IF( iom_use('diaharm_fast_eg_ssh_ts') ) THEN
167	! IF (lwp) CALL iom_put( 'diaharm_fast_eg_ssh_ts', sshn(5,5) )
168	!ENDIF
169
170	IF ( ln_diaharm_fast .and. ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
171
172	! this bit done every time step
173	nhm=2*nb_ana+1
174	c(1) = 1.0
175
176	sec2start_old = nint( (fjulday-fjulday_startharm)*86400._wp )
177	sec2start = nsec_day - NINT(0.5_wp * rdt)
178
179	sec2start = adatrj * 86400._wp
180
181
182	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,) 'diaharm_fast: sec2start = ',nint( (fjulday-fjulday_startharm)86400._wp ),nsec_day - NINT(0.5_wp * rdt),adatrj * 86400._wp
183
184	IF( iom_use('tide_t') ) CALL iom_put( 'tide_t', sec2start )
185
186
187
188
189
190	!IF(lwp) WRITE(numout,*) "ztime NEW", kt, sec2start, fjulday_startharm
191
192	DO jh=1,nb_ana
193	c(2jh ) = anaf(jh)cos( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
194	c(2jh+1) = anaf(jh)sin( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
195
196	c(2jh ) = anaf(jh)cos( sec2start*om_tide(jh) + anau(jh) )
197	c(2jh+1) = anaf(jh)sin( sec2start*om_tide(jh) + anau(jh) )
198
199
200
201	!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)
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, &
370	& ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d, &
371	& tname,ln_diaharm_multiyear,nn_diaharm_multiyear,ln_diaharm_update_nodal_daily,ln_diaharm_postproc_vel, ln_diaharm_verbose
372	!!----------------------------------------------------------------------
373	!JT
374	ln_diaharm_fast = .FALSE.
375	ln_diaharm_multiyear = .FALSE.
376	nn_diaharm_multiyear = 20
377	!JT
378	lk_diaharm_2D = .TRUE. ! to run 2d
379	lk_diaharm_3D = .TRUE. ! to run 3d
380
381	ln_diaharm_store = .TRUE.
382
383	ln_diaharm_postproc_vel = .FALSE.
384
385	IF(lwp) WRITE(numout,*)
386	IF(lwp) WRITE(numout,*) 'harm_init : initialization of harmonic analysis of tides'
387	IF(lwp) WRITE(numout,*) '~~~~~~~~~'
388
389	! GET NAMELIST DETAILS
390	REWIND( numnam_ref ) ! Namelist nam_diaharm_fast in reference namelist : Tidal harmonic analysis
391	READ ( numnam_ref, nam_diaharm_fast, IOSTAT = ios, ERR = 901)
392	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in reference namelist' )
393
394	REWIND( numnam_cfg ) ! Namelist nam_diaharm_fast in configuration namelist : Tidal harmonic analysis
395	READ ( numnam_cfg, nam_diaharm_fast, IOSTAT = ios, ERR = 902 )
396	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in configuration namelist' )
397	IF(lwm) WRITE ( numond, nam_diaharm_fast )
398
399
400	!
401	IF(lwp) THEN
402	WRITE(numout,*) 'Tidal diagnostics_fast'
403	WRITE(numout,*) ' Fast Harmonic analysis ?: ln_diaharm_fast= ', ln_diaharm_fast
404	WRITE(numout,*) ' Store output in restart?: ln_diaharm_store= ', ln_diaharm_store
405	WRITE(numout,*) ' Compute analysis?: ln_diaharm_compute= ', ln_diaharm_compute
406	WRITE(numout,*) ' Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
407	WRITE(numout,*) ' SSH harmonic analysis: ln_ana_ssh = ', ln_ana_ssh
408	WRITE(numout,*) ' Barotropic Velocities harmonic analysis: ln_ana_uvbar = ', ln_ana_uvbar
409	WRITE(numout,*) ' Bed Friction for harmonic analysis (not implemented): ln_ana_bfric = ', ln_ana_bfric
410	WRITE(numout,*) ' Density harmonic analysis: ln_ana_rho = ', ln_ana_rho
411	WRITE(numout,*) ' 3D velocities harmonic analysis: ln_ana_uv3d = ', ln_ana_uv3d
412	WRITE(numout,*) ' Vertical Velocities harmonic analysis: ln_ana_w3d = ', ln_ana_w3d
413	WRITE(numout,*) ' Names of harmonics: tname = ', tname
414	WRITE(numout,*) ' Max number of harmonics: jpmax_harmo = ', jpmax_harmo
415	WRITE(numout,*) ' Number of Harmonics: nb_harmo = ', nb_harmo
416	WRITE(numout,*) ' Multi-year harmonic analysis: ln_diaharm_multiyear = ', ln_diaharm_multiyear
417	WRITE(numout,*) ' Multi-year harmonic analysis - number of years: nn_diaharm_multiyear = ', nn_diaharm_multiyear
418	WRITE(numout,*) ' Multi-year harmonic analysis - number of years: ln_diaharm_update_nodal_daily = ', ln_diaharm_update_nodal_daily
419	WRITE(numout,*) ' Number of Harmonics: nyear, nmonth = ', nyear, nmonth
420	WRITE(numout,*) ' Post-process velocity stats: ln_diaharm_postproc_vel = ',ln_diaharm_postproc_vel
421	WRITE(numout,*) ' Verbose: ln_diaharm_verbose = ', ln_diaharm_verbose
422
423	ENDIF
424	! JT
425
426
427	IF ( ln_diaharm_multiyear ) THEN
428	ln_diaharm_store = .True.
429	IF(lwp) WRITE(numout,*) ' Multi-year harmonic analysis ', nyear,nn_diaharm_multiyear,nmonth
430	IF ((mod(nyear,nn_diaharm_multiyear) == 0) .AND. ( nmonth == 1 )) THEN ! Jan, year = 1980,2000,2020,2040, restart tidal calculation
431	ln_diaharm_read_restart = .FALSE.
432	IF(lwp) WRITE(numout,*) ' Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
433	ELSE
434	ln_diaharm_read_restart = .TRUE.
435	IF(lwp) WRITE(numout,*) ' Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
436	ENDIF
437
438
439
440	IF(lwp) WRITE(numout,*) ' Multi-year harmonic analysis ', nyear,nn_diaharm_multiyear,nmonth
441	IF ((mod(nyear,nn_diaharm_multiyear) == (nn_diaharm_multiyear - 1)) .AND. ( nmonth == 12 )) THEN ! Dec year = 1999,2019,2039,2040, restart tidal calculation
442	ln_diaharm_compute = .TRUE.
443	IF(lwp) WRITE(numout,*) ' Compute analysis?: ln_diaharm_compute= ', ln_diaharm_compute
444	ELSE
445	ln_diaharm_compute = .FALSE.
446	IF(lwp) WRITE(numout,*) ' Compute analysis?: ln_diaharm_compute= ', ln_diaharm_compute
447	ENDIF
448
449	ENDIF
450	IF ( kt < 10 ) THEN
451	ln_diaharm_read_restart = .FALSE.
452	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) ' kt = ',kt
453	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) ' kt < 10, so setting ln_diaharm_read_restart to .FALSE.'
454	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) ' Read in restart? : ln_diaharm_read_restart = ', ln_diaharm_read_restart
455	ENDIF
456
457	! JT
458
459
460	! GET NUMBER OF HARMONIC TO ANALYSE - from diaharm.F90
461	nb_ana = 0
462	DO jk=1,jpmax_harmo
463	DO ji=1,nb_harmo
464	IF(TRIM(tname(jk)) == TRIM(Wave( ntide(ji) )%cname_tide) ) THEN
465	nb_ana=nb_ana+1
466	ENDIF
467	END DO
468	END DO
469	!
470	IF(lwp) THEN
471	WRITE(numout,*) ' Namelist nam_diaharm_fast'
472	WRITE(numout,*) ' nb_ana = ', nb_ana
473	CALL flush(numout)
474	ENDIF
475	!
476	IF (nb_ana > nharm_max) THEN
477	IF(lwp) WRITE(numout,*) ' E R R O R harm_ana : nb_ana must be lower than nharm_max, stop'
478	IF(lwp) WRITE(numout,*) ' nharm_max = ', nharm_max
479	nstop = nstop + 1
480	ENDIF
481
482
483	ALLOCATE(ntide_all(nb_ana))
484	ALLOCATE(ntide_sub(nb_ana))
485
486	DO jk=1,nb_ana
487	DO ji=1,nb_harmo
488	IF (TRIM(tname(jk)) .eq. Wave( ntide(ji) )%cname_tide ) THEN
489	ntide_sub(jk) = ji
490	ntide_all(jk) = ntide(ji)
491	EXIT
492	END IF
493	END DO
494	END DO
495
496	ALLOCATE( anau(nb_ana) )
497	ALLOCATE( anav(nb_ana) )
498	ALLOCATE( anaf(nb_ana) )
499
500
501	IF( ln_diaharm_fast ) THEN
502
503	! SEARCH HOW MANY VARIABLES 2D AND 3D TO PROCESS
504	nvar_2d = 0; nvar_3d = 0
505	IF ( ln_ana_ssh ) nvar_2d = nvar_2d + 1 ! analysis elevation
506	IF ( ln_ana_uvbar ) nvar_2d = nvar_2d + 2 ! analysis depth-averaged velocity
507	IF ( ln_ana_bfric ) nvar_2d = nvar_2d + 2 ! analysis bottom friction
508
509	IF ( ln_ana_rho ) nvar_3d = nvar_3d + 1 ! analysis density
510	IF ( ln_ana_uv3d ) nvar_3d = nvar_3d + 2 ! analysis 3D horizontal velocities
511	IF ( ln_ana_w3d ) nvar_3d = nvar_3d + 1 ! analysis 3D vertical velocity
512
513	! CHECK IF SOMETHING TO RUN
514	IF ( nvar_2d .eq. 0 ) lk_diaharm_2D = .FALSE. ! no 2d to run
515	IF ( nvar_3d .eq. 0 ) lk_diaharm_3D = .FALSE. ! no 3d to run
516	! IF ( nvar_2d .gt. 0 .and. nvar_3d .gt. 0 ) lk_diaharm_fast = .FALSE.
517	! IF ( .NOT. ln_diaharm_store ) lk_diaharm_fast = .FALSE.
518
519	IF ( ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
520
521	! DO ALLOCATIONS
522	IF ( lk_diaharm_2D ) THEN
523	ALLOCATE( g_cumul_var2D(nb_ana*2+1,jpi,jpj, nvar_2d) )
524	ALLOCATE( g_cosamp2D( 0:nb_ana*2+1,jpi,jpj, nvar_2d) )
525	ALLOCATE( g_sinamp2D( 0:nb_ana*2+1,jpi,jpj, nvar_2d) )
526	ALLOCATE( g_out2D (jpi,jpj) )
527	ALLOCATE( h_out2D (jpi,jpj) )
528	ALLOCATE( m_posi_2d( nvar_2d ) ); m_posi_2d(:)=0
529	ENDIF
530
531	IF ( lk_diaharm_3D ) THEN
532	ALLOCATE( g_cumul_var3D(nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
533	ALLOCATE( g_cosamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
534	ALLOCATE( g_sinamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
535	ALLOCATE( g_out3D (jpi,jpj,jpk) )
536	ALLOCATE( h_out3D (jpi,jpj,jpk) )
537	ALLOCATE( m_posi_3d( nvar_3d ) ); m_posi_3d(:)=0
538	ENDIF
539
540	ALLOCATE( cc(nb_ana2+1,nb_ana2+1) )
541	ALLOCATE( a (nb_ana2+1,nb_ana2+1) )
542	ALLOCATE( bzz(nb_ana*2+1) )
543	ALLOCATE( x (nb_ana*2+1) )
544	ALLOCATE( c (nb_ana*2+1) )
545	! END ALLOCATE
546
547	! STORE INDEX OF WHAT TO PRODUCE DEPENDING ON ACTIVATED LOGICAL
548	! MAKES THINGS EASIER AND FASTER LATER
549	! !!! UGLY !!!
550	jh = 1; k2d = 0;
551	IF ( ln_ana_ssh ) THEN
552	k2d = k2d + 1; m_posi_2d(k2d) = jh
553	IF(lwp) WRITE(numout,*) " - ssh harmonic analysis activated (ln_ana_ssh)"
554	ENDIF
555	jh = jh + 1
556	IF ( ln_ana_uvbar ) THEN
557	k2d = k2d + 1; m_posi_2d(k2d) = jh
558	jh = jh + 1
559	k2d = k2d + 1; m_posi_2d(k2d) = jh
560	IF(lwp) WRITE(numout,*) " - barotropic currents harmonic analysis activated (ln_ana_uvbar)"
561	ELSE
562	jh = jh + 1
563	ENDIF
564	jh = jh + 1
565	IF ( ln_ana_bfric ) THEN
566	k2d = k2d + 1; m_posi_2d(k2d) = jh
567	jh = jh + 1;
568	k2d = k2d + 1; m_posi_2d(k2d) = jh
569	IF(lwp) WRITE(numout,*) " - bottom friction harmonic analysis activated (ln_ana_vbfr)"
570	ELSE
571	jh = jh + 1
572	ENDIF
573
574	! and for 3D
575	jh = 1; k3d = 0;
576	IF ( ln_ana_rho ) THEN
577	k3d = k3d + 1; m_posi_3d(k3d) = jh
578	IF(lwp) WRITE(numout,*) " - 3D density harmonic analysis activated (ln_ana_rho)"
579	ENDIF
580	jh = jh + 1
581	IF ( ln_ana_uv3d ) THEN
582	k3d = k3d + 1; m_posi_3d(k3d) = jh
583	jh = jh + 1
584	k3d = k3d + 1; m_posi_3d(k3d) = jh
585	IF(lwp) WRITE(numout,*) " - 3D horizontal currents harmonic analysis activated (ln_ana_uv3d)"
586	ELSE
587	jh = jh + 1
588	ENDIF
589	jh = jh + 1
590	IF ( ln_ana_w3d ) THEN
591	k3d = k3d + 1; m_posi_3d(k3d) = jh
592	IF(lwp) WRITE(numout,*) " - 3D vertical currents harmonic analysis activated (ln_ana_w3d)"
593	ENDIF
594
595	! SELECT AND STORE FREQUENCIES
596	IF(lwp) WRITE(numout,*) 'Analysed frequency : ',nb_ana ,'Frequency '
597	DO jh=1,nb_ana
598	om_tide(jh) = omega_tide( ntide_sub(jh) )
599	IF(lwp) WRITE(numout,) ' - ',tname(jh),' ',om_tide(jh), (2rpi/3600.)/om_tide(jh),"hr"
600	ENDDO
601
602	! READ RESTART IF
603	IF ( ln_diaharm_read_restart ) THEN
604	IF (lwp) WRITE(numout,*) "Reading previous harmonic data from previous run. kt = ",kt
605	! Need to read in bssh bz, cc anau anav and anaf
606	call harm_rst_read ! This reads in from the previous day
607	! Currrently the data in in assci format
608	ELSE
609
610	IF (lwp) WRITE(numout,*) "Starting harmonic analysis from Fresh. kt = ",kt
611
612	IF ( lk_diaharm_2D ) g_cumul_var2D(:,:,:,: ) = 0.0_wp
613	IF ( lk_diaharm_3D ) g_cumul_var3D(:,:,:,:,:) = 0.0_wp
614	cc = 0.0_wp
615	a (:,:) = 0.0_wp ! NB
616	bzz (:) = 0.0_wp
617	x (:) = 0.0_wp
618	c (:) = 0.0_wp
619	anau (:) = 0.0_wp
620	anav (:) = 0.0_wp
621	anaf (:) = 0.0_wp
622
623	DO jh = 1, nb_ana
624	anau(jh) = utide ( ntide_sub(jh) )
625	anav(jh) = v0tide( ntide_sub(jh) )
626	anaf(jh) = ftide ( ntide_sub(jh) )
627	END DO
628
629	fjulday_startharm=fjulday !Set this at very start and store
630	!JT this is a mistake - only works on daily cycles, should use fjulnsec_dayday
631
632	IF (lwp) THEN
633	WRITE(numout,*) '--------------------------'
634	WRITE(numout,*) ' - Output anaf for check'
635	WRITE(numout,*) 'ANA F', anaf
636	WRITE(numout,*) 'ANA U', anau
637	WRITE(numout,*) 'ANA V', anav
638	WRITE(numout,*) 'fjulday',fjulday
639	WRITE(numout,*) 'fjulday_startharm',fjulday_startharm
640	WRITE(numout,*) 'nsec_day',nsec_day
641	WRITE(numout,*) 'kt',kt
642	WRITE(numout,*) '--------------------------'
643	ENDIF
644
645	ENDIF
646
647	ELSE
648
649	IF (lwp) WRITE(numout,*) "No variable setup for harmonic analysis"
650
651	ENDIF
652	ENDIF
653
654	END SUBROUTINE harm_ana_init
655	!
656	SUBROUTINE gelim (a,b,x,n)
657	!!----------------------------------------------------------------------
658	!! * ROUTINE harm_ana *
659	!!
660	!! ** Purpose : Guassian elimination
661	!!
662	!!
663	!! ** Action : - first action (share memory array/varible modified
664	!! in this routine
665	!! - second action .....
666	!! - .....
667	!!
668	!! References :
669	!! Give references if exist otherwise suppress these lines
670	!!
671	!! History :
672	implicit none
673	!
674	integer :: n
675	REAL(WP) :: b(nb_ana2+1), a(nb_ana2+1,nb_ana*2+1)
676	REAL(WP) :: x(nb_ana*2+1)
677	INTEGER :: row,col,prow,pivrow,rrow
678	REAL(WP) :: atemp
679	REAL(WP) :: pivot
680	REAL(WP) :: m
681
682	do row=1,n-1
683	pivrow=row
684	pivot=a(row,n-row+1)
685	do prow=row+1,n
686	if (abs(a(prow,n-row+1)).gt.abs(pivot) ) then
687	pivot=a(prow,n-row+1)
688	pivrow=prow
689	endif
690	enddo
691	! swap row and prow
692	if ( pivrow .ne. row ) then
693	atemp=b(pivrow)
694	b(pivrow)=b(row)
695	b(row)=atemp
696	do col=1,n
697	atemp=a(pivrow,col)
698	a(pivrow,col)=a(row,col)
699	a(row,col)=atemp
700	enddo
701	endif
702
703	do rrow=row+1,n
704	if (a(row,row).ne.0) then
705
706	m=-a(rrow,n-row+1)/a(row,n-row+1)
707	do col=1,n
708	a(rrow,col)=m*a(row,col)+a(rrow,col)
709	enddo
710	b(rrow)=m*b(row)+b(rrow)
711	endif
712	enddo
713	enddo
714	! back substitution now
715
716	x(1)=b(n)/a(n,1)
717	do row=n-1,1,-1
718	x(n-row+1)=b(row)
719	do col=1,(n-row)
720	x(n-row+1)=(x(n-row+1)-a(row,col)*x(col))
721	enddo
722
723	x(n-row+1)=(x(n-row+1)/a(row,(n-row)+1))
724	enddo
725
726	return
727	END SUBROUTINE gelim
728
729	SUBROUTINE harm_ana_out
730	!!----------------------------------------------------------------------
731	!! * ROUTINE harm_ana_init *
732	!!
733	!! ** Purpose : initialization of ....
734	!!
735	!! ** Method : blah blah blah ...
736	!!
737	!! ** input : Namlist namexa
738	!!
739	!! ** Action : ...
740	!!
741	!! history :
742	!! 9.0 ! 03-08 (Autor Names) Original code
743	!!----------------------------------------------------------------------
744	USE dianam ! build name of file (routine)
745
746	!! * local declarations
747	INTEGER :: ji, jj, jk, jgrid, jh ! dummy loop indices
748	! INTEGER :: nh_T
749	! INTEGER :: nid_harm
750	! CHARACTER (len=40) :: clhstnamt, clop1, clop2 ! temporary names
751	! CHARACTER (len=40) :: clhstnamu, clhstnamv ! temporary names
752	CHARACTER (len=40) :: suffix
753	CHARACTER (len=40) :: tmp_name
754	! REAL(wp) :: zsto1, zsto2, zout, zmax, zjulian, zdt, zmdi ! temporary scalars
755
756	REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: amp_u2d,phi_u2d, amp_v2d,phi_v2d ! arrays for output
757	REAL(wp), ALLOCATABLE,DIMENSION(:,:,:,:) :: amp_u3d,phi_u3d, amp_v3d,phi_v3d ! arrays for output
758
759	REAL(wp) :: tmp_u_amp ,tmp_v_amp ,tmp_u_phi ,tmp_v_phi
760	REAL(wp) :: a_u, b_u, a_v, b_v, twodelta, delta, alpha2, alpha, qmin, qmax, ecc,thetamax, thetamin
761	REAL(wp) :: Qc, Qac, gc,gac, Phi_Ua, dir_Ua, polarity
762	REAL(wp) :: tmpreal
763
764	REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: tmp_u_amp_2d_mat,tmp_v_amp_2d_mat,tmp_u_phi_2d_mat,tmp_v_phi_2d_mat
765	REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: a_u_2d_mat,b_u_2d_mat,a_v_2d_mat,b_v_2d_mat
766	REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: qmax_2d_mat,qmin_2d_mat,ecc_2d_mat
767	REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: thetamax_2d_mat,thetamin_2d_mat,Qc_2d_mat,Qac_2d_mat
768	REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: gc_2d_mat,gac_2d_mat,Phi_Ua_2d_mat,dir_Ua_2d_mat
769	REAL(wp), ALLOCATABLE,DIMENSION(:,:) :: polarity_2d_mat
770
771	REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: tmp_u_amp_3d_mat,tmp_v_amp_3d_mat,tmp_u_phi_3d_mat,tmp_v_phi_3d_mat
772	REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: a_u_3d_mat,b_u_3d_mat,a_v_3d_mat,b_v_3d_mat
773	REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: qmax_3d_mat,qmin_3d_mat,ecc_3d_mat
774	REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: thetamax_3d_mat,thetamin_3d_mat,Qc_3d_mat,Qac_3d_mat
775	REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: gc_3d_mat,gac_3d_mat,Phi_Ua_3d_mat,dir_Ua_3d_mat
776	REAL(wp), ALLOCATABLE,DIMENSION(:,:,:) :: polarity_3d_mat
777
778
779
780
781	IF (ln_diaharm_postproc_vel) THEN
782	IF (ln_ana_uvbar) THEN
783	ALLOCATE( amp_u2d(nb_ana,jpi,jpj), amp_v2d(nb_ana,jpi,jpj), phi_u2d(nb_ana,jpi,jpj), phi_v2d(nb_ana,jpi,jpj) )
784
785	ALLOCATE(tmp_u_amp_2d_mat(jpi,jpj),tmp_v_amp_2d_mat(jpi,jpj),tmp_u_phi_2d_mat(jpi,jpj),tmp_v_phi_2d_mat(jpi,jpj))
786	ALLOCATE(a_u_2d_mat(jpi,jpj),b_u_2d_mat(jpi,jpj),a_v_2d_mat(jpi,jpj),b_v_2d_mat(jpi,jpj))
787	ALLOCATE(qmax_2d_mat(jpi,jpj),qmin_2d_mat(jpi,jpj),ecc_2d_mat(jpi,jpj))
788	ALLOCATE(thetamax_2d_mat(jpi,jpj),thetamin_2d_mat(jpi,jpj),Qc_2d_mat(jpi,jpj),Qac_2d_mat(jpi,jpj))
789	ALLOCATE(gc_2d_mat(jpi,jpj),gac_2d_mat(jpi,jpj),Phi_Ua_2d_mat(jpi,jpj),dir_Ua_2d_mat(jpi,jpj))
790	ALLOCATE(polarity_2d_mat(jpi,jpj))
791
792	ENDIF
793
794
795	IF (ln_ana_uv3d) THEN
796	ALLOCATE( amp_u3d(nb_ana,jpi,jpj,jpk), amp_v3d(nb_ana,jpi,jpj,jpk), phi_u3d(nb_ana,jpi,jpj,jpk), phi_v3d(nb_ana,jpi,jpj,jpk) )
797
798	ALLOCATE(tmp_u_amp_3d_mat(jpi,jpj,jpk),tmp_v_amp_3d_mat(jpi,jpj,jpk),tmp_u_phi_3d_mat(jpi,jpj,jpk),tmp_v_phi_3d_mat(jpi,jpj,jpk))
799	ALLOCATE(a_u_3d_mat(jpi,jpj,jpk),b_u_3d_mat(jpi,jpj,jpk),a_v_3d_mat(jpi,jpj,jpk),b_v_3d_mat(jpi,jpj,jpk))
800	ALLOCATE(qmax_3d_mat(jpi,jpj,jpk),qmin_3d_mat(jpi,jpj,jpk),ecc_3d_mat(jpi,jpj,jpk))
801	ALLOCATE(thetamax_3d_mat(jpi,jpj,jpk),thetamin_3d_mat(jpi,jpj,jpk),Qc_3d_mat(jpi,jpj,jpk),Qac_3d_mat(jpi,jpj,jpk))
802	ALLOCATE(gc_3d_mat(jpi,jpj,jpk),gac_3d_mat(jpi,jpj,jpk),Phi_Ua_3d_mat(jpi,jpj,jpk),dir_Ua_3d_mat(jpi,jpj,jpk))
803	ALLOCATE(polarity_3d_mat(jpi,jpj,jpk))
804
805	ENDIF
806	ENDIF
807
808	do jgrid=1,nvar_2d
809	do jh=1,nb_ana
810	h_out2D = 0.0
811	g_out2D = 0.0
812	do jj=1,nlcj
813	do ji=1,nlci
814	cca=g_cosamp2D(jh,ji,jj,jgrid)
815	ssa=g_sinamp2D(jh,ji,jj,jgrid)
816	h_out2D(ji,jj)=sqrt(cca2+ssa2)
817	IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
818	g_out2D(ji,jj)= 0.0_wp
819	ELSE
820	g_out2D(ji,jj)=(180.0/rpi)*atan2(ssa,cca)
821	ENDIF
822	IF (h_out2D(ji,jj).ne.0) THEN
823	h_out2D(ji,jj)=h_out2D(ji,jj)/anaf(jh)
824	ENDIF
825	IF (g_out2D(ji,jj).ne.0) THEN !Correct and take modulus
826	!JT
827	!JT g_out2D(ji,jj) = g_out2D(ji,jj) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
828	!JT
829	g_out2D(ji,jj) = g_out2D(ji,jj) + MOD( (anau(jh))/rad , 360.0)
830	if (g_out2D(ji,jj).gt.360.0) then
831	g_out2D(ji,jj)=g_out2D(ji,jj)-360.0
832	else if (g_out2D(ji,jj).lt.0.0) then
833	g_out2D(ji,jj)=g_out2D(ji,jj)+360.0
834	endif
835	ENDIF
836	enddo
837	enddo
838	!
839	! NETCDF OUTPUT
840	suffix = TRIM( m_varName2d( m_posi_2d(jgrid) ) )
841	IF(lwp) WRITE(numout,*) "diaharm_fast", suffix
842
843	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix)
844	IF( iom_use(TRIM(tmp_name)) ) THEN
845	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out2D)
846	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,) "diaharm_fast names", tmp_name,tname(jh),' ',om_tide(jh), (2rpi/3600.)/om_tide(jh),"hr"
847	CALL iom_put( TRIM(tmp_name), h_out2D(:,:) )
848	ELSE
849	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
850	ENDIF
851
852	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix)
853	IF( iom_use(TRIM(tmp_name)) ) THEN
854	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out2D)
855	CALL iom_put( TRIM(tmp_name), g_out2D(:,:) )
856	ELSE
857	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
858	ENDIF
859
860
861
862	IF (ln_diaharm_postproc_vel .AND. ln_ana_uvbar) THEN
863
864	!IF (m_posi_2d(jgrid) == 2) THEN
865	IF (TRIM(suffix) == TRIM('u2d')) THEN
866	if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' u2d '//TRIM(suffix)
867	do jj=1,nlcj
868	do ji=1,nlci
869	if (ssumask(ji,jj) == 1) THEN
870	amp_u2d(jh,ji,jj) = h_out2D(ji,jj)
871	phi_u2d(jh,ji,jj) = rpi*g_out2D(ji,jj)/180.0
872	else
873	amp_u2d(jh,ji,jj) = 0.
874	phi_u2d(jh,ji,jj) = 0.
875	ENDIF
876	enddo
877	enddo
878	ENDIF
879
880	!IF (m_posi_2d(jgrid) == 3) THEN
881	IF (TRIM(suffix) == TRIM('v2d')) THEN
882	if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' v2d '//TRIM(suffix)
883	do jj=1,nlcj
884	do ji=1,nlci
885	if (ssvmask(ji,jj) == 1) THEN
886	amp_v2d(jh,ji,jj) = h_out2D(ji,jj)
887	phi_v2d(jh,ji,jj) = rpi*g_out2D(ji,jj)/180.0
888	else
889	amp_v2d(jh,ji,jj) = 0
890	phi_v2d(jh,ji,jj) = 0
891	ENDIF
892	enddo
893	enddo
894	ENDIF
895	ENDIF
896
897	CALL FLUSH(numout)
898
899
900	enddo
901
902	suffix = TRIM( m_varName2d( m_posi_2d(jgrid) ) )
903	tmp_name='TA_'//TRIM(suffix)//'_off'
904	IF( iom_use(TRIM(tmp_name)) ) THEN
905	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
906	CALL iom_put( TRIM(tmp_name), g_cosamp2D( 0,:,:,jgrid))
907	ELSE
908	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
909	ENDIF
910
911	CALL FLUSH(numout)
912
913	enddo
914	!
915	! DO THE SAME FOR 3D VARIABLES
916	!
917	do jgrid=1,nvar_3d
918	do jh=1,nb_ana
919	h_out3D = 0.0
920	g_out3D = 0.0
921	DO jk=1,jpkm1
922	do jj=1,nlcj
923	do ji=1,nlci
924	cca=g_cosamp3D(jh,ji,jj,jk,jgrid)
925	ssa=g_sinamp3D(jh,ji,jj,jk,jgrid)
926	h_out3D(ji,jj,jk)=sqrt(cca2+ssa2)
927	IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
928	g_out3D(ji,jj,jk) = 0.0_wp
929	ELSE
930	g_out3D(ji,jj,jk) = (180.0/rpi)*atan2(ssa,cca)
931	ENDIF
932	IF (h_out3D(ji,jj,jk).ne.0) THEN
933	h_out3D(ji,jj,jk) = h_out3D(ji,jj,jk)/anaf(jh)
934	ENDIF
935	IF (g_out3D(ji,jj,jk).ne.0) THEN !Correct and take modulus
936	!JT
937	!JT g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
938	!JT
939	g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk) + MOD( (anau(jh))/rad , 360.0)
940	if (g_out3D(ji,jj,jk).gt.360.0) then
941	g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)-360.0
942	else if (g_out3D(ji,jj,jk).lt.0.0) then
943	g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)+360.0
944	endif
945	ENDIF
946	enddo ! ji
947	enddo ! jj
948	ENDDO ! jk
949	!
950	! NETCDF OUTPUT
951	suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) )
952	IF(lwp) WRITE(numout,*) "diaharm_fast", suffix
953
954	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix)
955	IF( iom_use(TRIM(tmp_name)) ) THEN
956	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(h_out3D)
957	CALL iom_put( TRIM(tmp_name), h_out3D(:,:,:) )
958	ELSE
959	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
960	ENDIF
961
962	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix)
963	IF( iom_use(TRIM(tmp_name)) ) THEN
964	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(g_out3D)
965	CALL iom_put(tmp_name, g_out3D(:,:,:) )
966	ELSE
967	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
968	ENDIF
969
970
971
972	IF (ln_diaharm_postproc_vel .AND. ln_ana_uv3d) THEN
973
974	!IF (m_posi_2d(jgrid) == 2) THEN
975	IF (TRIM(suffix) == TRIM('u3d')) THEN
976	if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' u3d '//TRIM(suffix)
977	DO jk=1,jpkm1
978	do jj=1,nlcj
979	do ji=1,nlci
980	if (umask(ji,jj,jk) == 1) THEN
981	amp_u3d(jh,ji,jj,jk) = h_out3D(ji,jj,jk)
982	phi_u3d(jh,ji,jj,jk) = rpi*g_out3D(ji,jj,jk)/180.0
983	else
984	amp_u3d(jh,ji,jj,jk) = 0.
985	phi_u3d(jh,ji,jj,jk) = 0.
986	ENDIF
987	enddo
988	enddo
989	enddo
990	ENDIF
991
992	!IF (m_posi_2d(jgrid) == 3) THEN
993	IF (TRIM(suffix) == TRIM('v3d')) THEN
994	if (lwp) WRITE(numout,*) "diaharm_fast ln_diaharm_postproc_vel: "//TRIM(Wave(ntide_all(jh))%cname_tide)//' v3d '//TRIM(suffix)
995	DO jk=1,jpkm1
996	do jj=1,nlcj
997	do ji=1,nlci
998	if (vmask(ji,jj,jk) == 1) THEN
999	amp_v3d(jh,ji,jj,jk) = h_out3D(ji,jj,jk)
1000	phi_v3d(jh,ji,jj,jk) = rpi*g_out3D(ji,jj,jk)/180.0
1001	else
1002	amp_v3d(jh,ji,jj,jk) = 0
1003	phi_v3d(jh,ji,jj,jk) = 0
1004	ENDIF
1005	enddo
1006	enddo
1007	enddo
1008	ENDIF
1009	ENDIF
1010
1011	CALL FLUSH(numout)
1012
1013
1014	enddo ! jh
1015
1016	suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) )
1017	tmp_name='TA_'//TRIM(suffix)//'_off'
1018	IF( iom_use(TRIM(tmp_name)) ) THEN
1019	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1020	CALL iom_put( TRIM(tmp_name), g_cosamp3D( 0,:,:,:,jgrid))
1021	ELSE
1022	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: not requested: ",TRIM(tmp_name)
1023	ENDIF
1024
1025	enddo ! jgrid
1026
1027	CALL FLUSH(numout)
1028
1029
1030	IF (ln_diaharm_postproc_vel ) THEN
1031	IF ( ln_ana_uvbar) THEN
1032	IF(lwp) WRITE(numout,*) "diaharm_fast: Postprocess barotropic velocity tidal parameters"
1033	CALL FLUSH(numout)
1034	DO jh=1,nb_ana
1035
1036
1037	tmp_u_amp_2d_mat(:,:) = 0.
1038	tmp_v_amp_2d_mat(:,:) = 0.
1039	tmp_u_phi_2d_mat(:,:) = 0.
1040	tmp_v_phi_2d_mat(:,:) = 0.
1041
1042	a_u_2d_mat(:,:) = 0.
1043	b_u_2d_mat(:,:) = 0.
1044	a_v_2d_mat(:,:) = 0.
1045	b_v_2d_mat(:,:) = 0.
1046
1047	qmax_2d_mat(:,:) = 0.
1048	qmin_2d_mat(:,:) = 0.
1049
1050	ecc_2d_mat(:,:) = 0.
1051	thetamax_2d_mat(:,:) =0.
1052	thetamin_2d_mat(:,:) = 0.
1053
1054	Qc_2d_mat(:,:) = 0.
1055	Qac_2d_mat(:,:) = 0.
1056	gc_2d_mat(:,:) = 0.
1057	gac_2d_mat(:,:) = 0.
1058
1059	Phi_Ua_2d_mat(:,:) = 0.
1060	dir_Ua_2d_mat(:,:) = 0.
1061	polarity_2d_mat(:,:) = 0.
1062
1063
1064	DO jj = 1, nlcj !- 1
1065	DO ji = 1, nlci ! - 1
1066
1067	IF ( ((ssumask(ji,jj) + ssumask(ji-1,jj)) > 0 ) .AND. ((ssvmask(ji,jj) + ssvmask(ji,jj-1)) > 0 ) ) THEN
1068
1069	if ( (ssumask(ji,jj) == 1) .AND. (ssumask(ji-1,jj) == 1)) then
1070
1071	tmp_u_amp = ((amp_u2d(jh,ji,jj)ssumask(ji,jj)) + (amp_u2d(jh,ji-1,jj)ssumask(ji-1,jj)))/(ssumask(ji,jj) + ssumask(ji-1,jj))
1072	!tmp_u_phi = ((phi_u2d(jh,ji,jj)ssumask(ji,jj)) + (phi_u2d(jh,ji-1,jj)ssumask(ji-1,jj)))/(ssumask(ji,jj) + ssumask(ji-1,jj))
1073	tmp_u_phi = atan2((sin(phi_u2d(jh,ji,jj)) + sin(phi_u2d(jh,ji-1,jj))),(cos(phi_u2d(jh,ji,jj)) + cos(phi_u2d(jh,ji-1,jj))))
1074	else if ( (ssumask(ji,jj) == 1) .AND. (ssumask(ji-1,jj) == 0)) then
1075	tmp_u_amp = (amp_u2d(jh,ji,jj)*ssumask(ji,jj))
1076	tmp_u_phi = (phi_u2d(jh,ji,jj)*ssumask(ji,jj))
1077	else if ( (ssumask(ji,jj) == 0) .AND. (ssumask(ji-1,jj) == 1)) then
1078	tmp_u_amp = (amp_u2d(jh,ji-1,jj)*ssumask(ji-1,jj))
1079	tmp_u_phi = (phi_u2d(jh,ji-1,jj)*ssumask(ji-1,jj))
1080	else
1081	cycle
1082	end if
1083
1084
1085	if ( (ssvmask(ji,jj) == 1) .AND. (ssvmask(ji,jj-1) == 1)) then
1086	tmp_v_amp = ((amp_v2d(jh,ji,jj)ssvmask(ji,jj)) + (amp_v2d(jh,ji,jj-1)ssvmask(ji,jj-1)))/(ssvmask(ji,jj) + ssvmask(ji,jj-1))
1087	!tmp_v_phi = ((phi_v2d(jh,ji,jj)ssvmask(ji,jj)) + (phi_v2d(jh,ji,jj-1)ssvmask(ji,jj-1)))/(ssvmask(ji,jj) + ssvmask(ji,jj-1))
1088	tmp_v_phi = atan2((sin(phi_v2d(jh,ji,jj)) + sin(phi_v2d(jh,ji,jj-1))),(cos(phi_v2d(jh,ji,jj)) + cos(phi_v2d(jh,ji,jj-1))))
1089	else if ( (ssvmask(ji,jj) == 1) .AND. (ssvmask(ji,jj-1) == 0)) then
1090	tmp_v_amp = (amp_v2d(jh,ji,jj)*ssvmask(ji,jj))
1091	tmp_v_phi = (phi_v2d(jh,ji,jj)*ssvmask(ji,jj))
1092	else if ( (ssvmask(ji,jj) == 0) .AND. (ssvmask(ji,jj-1) == 1)) then
1093	tmp_v_amp = (amp_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1))
1094	!tmp_v_phi = (phi_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1))
1095	tmp_v_phi = (phi_v2d(jh,ji,jj-1)*ssvmask(ji,jj-1))
1096	else
1097	cycle
1098	end if
1099
1100
1101	a_u = tmp_U_amp * cos(tmp_U_phi)
1102	b_u = tmp_U_amp * sin(tmp_U_phi)
1103	a_v = tmp_V_amp * cos(tmp_V_phi)
1104	b_v = tmp_V_amp * sin(tmp_V_phi)
1105
1106	twodelta = atan2( (tmp_V_amp*2 sin( 2(tmp_U_phi - tmp_V_phi) ) ) , ( tmp_U_amp2 + tmp_V_amp2 cos( 2*(tmp_U_phi - tmp_V_phi) ) ) )
1107	delta = twodelta/2.
1108
1109	!alpha2 = sqrt( tmp_U_amp4 + tmp_V_amp4 + 2tmp_U_amp2tmp_V_amp*2cos(2*(tmp_U_phi - tmp_V_phi)) )
1110	tmpreal = tmp_U_amp4 + tmp_V_amp4 + 2tmp_U_amp2tmp_V_amp*2cos(2*(tmp_U_phi - tmp_V_phi))
1111	if (tmpreal < 0) tmpreal = 0 !CYCLE
1112
1113	!alpha2 = sqrt( tmp_U_amp4 + tmp_V_amp4 + 2tmp_U_amp2tmp_V_amp*2cos(2*(tmp_U_phi - tmp_V_phi)) )
1114	alpha2 = sqrt( tmpreal )
1115	if (alpha2 < 0) alpha2 = 0 !CYCLE
1116	alpha= sqrt( alpha2 )
1117
1118
1119	!major and minor axis of the ellipse
1120	qmax = sqrt( (tmp_U_amp2 + tmp_V_amp2 + alpha**2)/2 )
1121
1122	tmpreal = (tmp_U_amp2 + tmp_V_amp2 - alpha**2)/2
1123	if (tmpreal < 0) tmpreal = 0 !CYCLE
1124	!qmin = sqrt( (tmp_U_amp2 + tmp_V_amp2 - alpha**2)/2 ) ! but always positive.
1125	qmin = sqrt( tmpreal ) ! but always positive.
1126
1127	!eccentricity of ellipse
1128	tmpreal = (qmax + qmin)
1129	if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE
1130	!ecc = (qmax - qmin)/(qmax + qmin)
1131	ecc = (qmax - qmin)/(tmpreal)
1132
1133	! Angle of major and minor ellipse
1134	thetamax = atan2(( tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta) ) , ( tmp_U_amp * cos( delta) ) )
1135	thetamin = thetamax + rpi/2.
1136
1137
1138
1139	! Rotary current components: Pugh A3.10
1140	! Clockwise (c) and anticlockwise (ac) rotating rotate_wind_vectors
1141	! so Qc = clockwise = anticyclonic = negative
1142	! and Qac = anticlockwise = cyclonic = negative
1143
1144	tmpreal = tmp_U_amp2 + tmp_V_amp2 - (2tmp_U_amptmp_V_amp*sin( tmp_V_phi - tmp_U_phi))
1145	if (tmpreal < 0) tmpreal = 0 !CYCLE
1146	!Qc = 0.5sqrt( tmp_U_amp2 + tmp_V_amp2 - (2tmp_U_amptmp_V_ampsin( tmp_V_phi - tmp_U_phi)) )
1147	Qc = 0.5*sqrt( tmpreal )
1148
1149	tmpreal = tmp_U_amp2 + tmp_V_amp2 + (2tmp_U_amptmp_V_amp*sin( tmp_V_phi - tmp_U_phi))
1150	if (tmpreal < 0) tmpreal = 0 !CYCLE
1151	!Qac = 0.5sqrt( tmp_U_amp2 + tmp_V_amp2 + (2tmp_U_amptmp_V_ampsin( tmp_V_phi - tmp_U_phi)) )
1152	Qac = 0.5*sqrt( tmpreal )
1153
1154
1155	gc = atan2( ( ( tmp_U_ampsin( tmp_U_phi ) ) + (tmp_V_ampcos( tmp_V_phi) ) ) , &
1156	& ( (tmp_U_ampcos( tmp_U_phi )) - (tmp_V_ampsin( tmp_V_phi )) ) )
1157	gac = atan2( ( ( -tmp_U_ampsin( tmp_U_phi ) ) + (tmp_V_ampcos( tmp_V_phi) ) ) , &
1158	& ( (tmp_U_ampcos( tmp_U_phi )) + (tmp_V_ampsin( tmp_V_phi )) ) )
1159
1160	!Pugh A3.2
1161	Phi_Ua = -0.5*(gac - gc)
1162	dir_Ua = 0.5*(gac + gc) ! positive from x axis
1163
1164	tmpreal = qmax
1165	!if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE
1166	if (tmpreal == 0) then
1167	polarity = 0
1168	else
1169	polarity = (Qac - Qc)/qmax
1170	endif
1171
1172	tmp_u_amp_2d_mat(ji,jj) = tmp_u_amp
1173	tmp_v_amp_2d_mat(ji,jj) = tmp_v_amp
1174	tmp_u_phi_2d_mat(ji,jj) = tmp_u_phi
1175	tmp_v_phi_2d_mat(ji,jj) = tmp_v_phi
1176
1177
1178	a_u_2d_mat(ji,jj) = a_u
1179	b_u_2d_mat(ji,jj) = b_u
1180	a_v_2d_mat(ji,jj) = a_v
1181	b_v_2d_mat(ji,jj) = b_v
1182
1183	qmax_2d_mat(ji,jj) = qmax
1184	qmin_2d_mat(ji,jj) = qmin
1185
1186	ecc_2d_mat(ji,jj) = ecc
1187	thetamax_2d_mat(ji,jj) = thetamax
1188	thetamin_2d_mat(ji,jj) = thetamin
1189
1190	Qc_2d_mat(ji,jj) = Qc
1191	Qac_2d_mat(ji,jj) = Qac
1192	gc_2d_mat(ji,jj) = gc
1193	gac_2d_mat(ji,jj) = gac
1194
1195	Phi_Ua_2d_mat(ji,jj) = Phi_Ua
1196	dir_Ua_2d_mat(ji,jj) = dir_Ua
1197	polarity_2d_mat(ji,jj) = polarity
1198
1199	ENDIF
1200	END DO
1201	END DO
1202
1203
1204	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_amp_t_uvbar'
1205	IF( iom_use(TRIM(tmp_name)) ) THEN
1206	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1207	CALL iom_put( TRIM(tmp_name), tmp_u_amp_2d_mat(:,:))
1208	ENDIF
1209	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_amp_t_uvbar'
1210	IF( iom_use(TRIM(tmp_name)) ) THEN
1211	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1212	CALL iom_put( TRIM(tmp_name), tmp_v_amp_2d_mat(:,:))
1213	ENDIF
1214	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_phi_t_uvbar'
1215	IF( iom_use(TRIM(tmp_name)) ) THEN
1216	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1217	CALL iom_put( TRIM(tmp_name), tmp_u_phi_2d_mat(:,:))
1218	ENDIF
1219	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_phi_t_uvbar'
1220	IF( iom_use(TRIM(tmp_name)) ) THEN
1221	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1222	CALL iom_put( TRIM(tmp_name), tmp_v_phi_2d_mat(:,:))
1223	ENDIF
1224
1225
1226
1227	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_u_uvbar'
1228	IF( iom_use(TRIM(tmp_name)) ) THEN
1229	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1230	CALL iom_put( TRIM(tmp_name), a_u_2d_mat(:,:))
1231	ENDIF
1232	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_v_uvbar'
1233	IF( iom_use(TRIM(tmp_name)) ) THEN
1234	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1235	CALL iom_put( TRIM(tmp_name), a_v_2d_mat(:,:))
1236	ENDIF
1237	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_u_uvbar'
1238	IF( iom_use(TRIM(tmp_name)) ) THEN
1239	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1240	CALL iom_put( TRIM(tmp_name), b_u_2d_mat(:,:))
1241	ENDIF
1242	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_v_uvbar'
1243	IF( iom_use(TRIM(tmp_name)) ) THEN
1244	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1245	CALL iom_put( TRIM(tmp_name), b_v_2d_mat(:,:))
1246	ENDIF
1247
1248	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmax_uvbar'
1249	IF( iom_use(TRIM(tmp_name)) ) THEN
1250	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1251	CALL iom_put( TRIM(tmp_name), qmax_2d_mat(:,:))
1252	ENDIF
1253	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmin_uvbar'
1254	IF( iom_use(TRIM(tmp_name)) ) THEN
1255	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1256	CALL iom_put( TRIM(tmp_name), qmin_2d_mat(:,:))
1257	ENDIF
1258
1259	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_ecc_uvbar'
1260	IF( iom_use(TRIM(tmp_name)) ) THEN
1261	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1262	CALL iom_put( TRIM(tmp_name), ecc_2d_mat(:,:))
1263	ENDIF
1264	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamax_uvbar'
1265	IF( iom_use(TRIM(tmp_name)) ) THEN
1266	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1267	CALL iom_put( TRIM(tmp_name), thetamax_2d_mat(:,:))
1268	ENDIF
1269	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamin_uvbar'
1270	IF( iom_use(TRIM(tmp_name)) ) THEN
1271	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1272	CALL iom_put( TRIM(tmp_name), thetamin_2d_mat(:,:))
1273	ENDIF
1274
1275	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qc_uvbar'
1276	IF( iom_use(TRIM(tmp_name)) ) THEN
1277	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1278	CALL iom_put( TRIM(tmp_name), Qc_2d_mat(:,:))
1279	ENDIF
1280	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qac_uvbar'
1281	IF( iom_use(TRIM(tmp_name)) ) THEN
1282	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1283	CALL iom_put( TRIM(tmp_name), Qac_2d_mat(:,:))
1284	ENDIF
1285	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gc_uvbar'
1286	IF( iom_use(TRIM(tmp_name)) ) THEN
1287	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1288	CALL iom_put( TRIM(tmp_name), gc_2d_mat(:,:))
1289	ENDIF
1290	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gac_uvbar'
1291	IF( iom_use(TRIM(tmp_name)) ) THEN
1292	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1293	CALL iom_put( TRIM(tmp_name), gac_2d_mat(:,:))
1294	ENDIF
1295
1296
1297	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Phi_Ua_uvbar'
1298	IF( iom_use(TRIM(tmp_name)) ) THEN
1299	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1300	CALL iom_put( TRIM(tmp_name), Phi_Ua_2d_mat(:,:))
1301	ENDIF
1302	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_dir_Ua_uvbar'
1303	IF( iom_use(TRIM(tmp_name)) ) THEN
1304	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1305	CALL iom_put( TRIM(tmp_name), dir_Ua_2d_mat(:,:))
1306	ENDIF
1307	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_polarity_uvbar'
1308	IF( iom_use(TRIM(tmp_name)) ) THEN
1309	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1310	CALL iom_put( TRIM(tmp_name), polarity_2d_mat(:,:))
1311	ENDIF
1312
1313	tmp_u_amp_2d_mat(:,:) = 0.
1314	tmp_v_amp_2d_mat(:,:) = 0.
1315	tmp_u_phi_2d_mat(:,:) = 0.
1316	tmp_v_phi_2d_mat(:,:) = 0.
1317
1318	a_u_2d_mat(:,:) = 0.
1319	b_u_2d_mat(:,:) = 0.
1320	a_v_2d_mat(:,:) = 0.
1321	b_v_2d_mat(:,:) = 0.
1322
1323	qmax_2d_mat(:,:) = 0.
1324	qmin_2d_mat(:,:) = 0.
1325
1326	ecc_2d_mat(:,:) = 0.
1327	thetamax_2d_mat(:,:) =0.
1328	thetamin_2d_mat(:,:) = 0.
1329
1330	Qc_2d_mat(:,:) = 0.
1331	Qac_2d_mat(:,:) = 0.
1332	gc_2d_mat(:,:) = 0.
1333	gac_2d_mat(:,:) = 0.
1334
1335	Phi_Ua_2d_mat(:,:) = 0.
1336	dir_Ua_2d_mat(:,:) = 0.
1337	polarity_2d_mat(:,:) = 0.
1338
1339
1340	END DO
1341	IF(lwp) WRITE(numout,*) "diaharm_fast: Finshed postprocessing 2d velocity tidal parameters"
1342	ENDIF
1343
1344	CALL FLUSH(numout)
1345
1346	IF (ln_ana_uv3d) THEN
1347	IF(lwp) WRITE(numout,*) "diaharm_fast: Postprocess baroclinic velocity tidal parameters"
1348	CALL FLUSH(numout)
1349	DO jh=1,nb_ana
1350
1351
1352	tmp_u_amp_3d_mat(:,:,:) = 0.
1353	tmp_v_amp_3d_mat(:,:,:) = 0.
1354	tmp_u_phi_3d_mat(:,:,:) = 0.
1355	tmp_v_phi_3d_mat(:,:,:) = 0.
1356
1357	a_u_3d_mat(:,:,:) = 0.
1358	b_u_3d_mat(:,:,:) = 0.
1359	a_v_3d_mat(:,:,:) = 0.
1360	b_v_3d_mat(:,:,:) = 0.
1361
1362	qmax_3d_mat(:,:,:) = 0.
1363	qmin_3d_mat(:,:,:) = 0.
1364
1365	ecc_3d_mat(:,:,:) = 0.
1366	thetamax_3d_mat(:,:,:) =0.
1367	thetamin_3d_mat(:,:,:) = 0.
1368
1369	Qc_3d_mat(:,:,:) = 0.
1370	Qac_3d_mat(:,:,:) = 0.
1371	gc_3d_mat(:,:,:) = 0.
1372	gac_3d_mat(:,:,:) = 0.
1373
1374	Phi_Ua_3d_mat(:,:,:) = 0.
1375	dir_Ua_3d_mat(:,:,:) = 0.
1376	polarity_3d_mat(:,:,:) = 0.
1377
1378
1379	DO jk=1,jpkm1
1380	!DO jj = 2, nlcj ! - 1
1381	! DO ji = 2, nlci ! - 1
1382
1383	! DO jj = 1, jpjm1 #works
1384	! DO ji = 1, jpim1
1385
1386	DO jj = 1, nlcj !- 1
1387	DO ji = 1, nlci ! - 1
1388
1389
1390	! do jj=2,nlcj
1391	! do ji=2,nlci
1392	!IF ((umask(ji,jj) + umask(ji-1,jj)) == 0 ) CYCLE
1393	!IF ((vmask(ji,jj) + vmask(ji,jj-1)) == 0 ) CYCLE
1394
1395	IF ( ((umask(ji,jj,jk) + umask(ji-1,jj,jk)) > 0 ) .AND. ((vmask(ji,jj,jk) + vmask(ji,jj-1,jk)) > 0 ) ) THEN
1396	!tmp_u_amp = ((amp_u3d(jh,ji,jj,jk)umask(ji,jj,jk)) + (amp_u3d(jh,ji-1,jj,jk)umask(ji-1,jj,jk)))/(umask(ji,jj,jk) + umask(ji-1,jj,jk))
1397	!tmp_v_amp = ((amp_v3d(jh,ji,jj,jk)vmask(ji,jj,jk)) + (amp_v3d(jh,ji,jj-1,jk)vmask(ji,jj-1,jk)))/(vmask(ji,jj,jk) + vmask(ji,jj-1,jk))
1398	! WORK ON THE WRAP AROUND
1399	!tmp_u_phi = ((phi_u3d(jh,ji,jj,jk)umask(ji,jj,jk)) + (phi_u3d(jh,ji-1,jj,jk)umask(ji-1,jj,jk)))/(umask(ji,jj,jk) + umask(ji-1,jj,jk))
1400	!tmp_v_phi = ((phi_v3d(jh,ji,jj,jk)vmask(ji,jj,jk)) + (phi_v3d(jh,ji,jj-1,jk)vmask(ji,jj-1,jk)))/(vmask(ji,jj,jk) + vmask(ji,jj-1,jk))
1401
1402	if ( (umask(ji,jj,jk) == 1) .AND. (umask(ji-1,jj,jk) == 1)) then
1403	tmp_u_amp = ((amp_u3d(jh,ji,jj,jk)umask(ji,jj,jk)) + (amp_u3d(jh,ji-1,jj,jk)umask(ji-1,jj,jk)))/(umask(ji,jj,jk) + umask(ji-1,jj,jk))
1404	!tmp_u_phi = ((phi_u3d(jh,ji,jj,jk)umask(ji,jj,jk)) + (phi_u3d(jh,ji-1,jj,jk)umask(ji-1,jj,jk)))/(umask(ji,jj,jk) + umask(ji-1,jj,jk))
1405	tmp_u_phi = atan2((sin(phi_u3d(jh,ji,jj,jk)) + sin(phi_u3d(jh,ji-1,jj,jk))),(cos(phi_u3d(jh,ji,jj,jk)) + cos(phi_u3d(jh,ji-1,jj,jk))))
1406	else if ( (umask(ji,jj,jk) == 1) .AND. (umask(ji-1,jj,jk) == 0)) then
1407	tmp_u_amp = (amp_u3d(jh,ji,jj,jk)*umask(ji,jj,jk))
1408	tmp_u_phi = (phi_u3d(jh,ji,jj,jk)*umask(ji,jj,jk))
1409	else if ( (umask(ji,jj,jk) == 0) .AND. (umask(ji-1,jj,jk) == 1)) then
1410	tmp_u_amp = (amp_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk))
1411	tmp_u_phi = (phi_u3d(jh,ji-1,jj,jk)*umask(ji-1,jj,jk))
1412	else
1413	cycle
1414	end if
1415
1416
1417	if ( (vmask(ji,jj,jk) == 1) .AND. (vmask(ji,jj-1,jk) == 1)) then
1418	tmp_v_amp = ((amp_v3d(jh,ji,jj,jk)vmask(ji,jj,jk)) + (amp_v3d(jh,ji,jj-1,jk)vmask(ji,jj-1,jk)))/(vmask(ji,jj,jk) + vmask(ji,jj-1,jk))
1419	!tmp_v_phi = ((phi_v3d(jh,ji,jj,jk)vmask(ji,jj,jk)) + (phi_v3d(jh,ji,jj-1,jk)vmask(ji,jj-1,jk)))/(vmask(ji,jj,jk) + vmask(ji,jj-1,jk))
1420	tmp_v_phi = atan2((sin(phi_v3d(jh,ji,jj,jk)) + sin(phi_v3d(jh,ji,jj-1,jk))),(cos(phi_v3d(jh,ji,jj,jk)) + cos(phi_v3d(jh,ji,jj-1,jk))))
1421	else if ( (vmask(ji,jj,jk) == 1) .AND. (vmask(ji,jj-1,jk) == 0)) then
1422	tmp_v_amp = (amp_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk))
1423	tmp_v_phi = (phi_v3d(jh,ji,jj,jk)*vmask(ji,jj,jk))
1424	else if ( (vmask(ji,jj,jk) == 0) .AND. (vmask(ji,jj-1,jk) == 1)) then
1425	tmp_v_amp = (amp_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk))
1426	tmp_v_phi = (phi_v3d(jh,ji,jj-1,jk)*vmask(ji,jj-1,jk))
1427	else
1428	cycle
1429	end if
1430
1431
1432
1433	a_u = tmp_U_amp * cos(tmp_U_phi)
1434	b_u = tmp_U_amp * sin(tmp_U_phi)
1435	a_v = tmp_V_amp * cos(tmp_V_phi)
1436	b_v = tmp_V_amp * sin(tmp_V_phi)
1437
1438	twodelta = atan2( (tmp_V_amp*2 sin( 2(tmp_U_phi - tmp_V_phi) ) ) , ( tmp_U_amp2 + tmp_V_amp2 cos( 2*(tmp_U_phi - tmp_V_phi) ) ) )
1439	delta = twodelta/2.
1440
1441	!alpha2 = sqrt( tmp_U_amp4 + tmp_V_amp4 + 2tmp_U_amp2tmp_V_amp*2cos(2*(tmp_U_phi - tmp_V_phi)) )
1442	tmpreal = tmp_U_amp4 + tmp_V_amp4 + 2tmp_U_amp2tmp_V_amp*2cos(2*(tmp_U_phi - tmp_V_phi))
1443	if (tmpreal < 0) tmpreal = 0 !CYCLE
1444
1445	!alpha2 = sqrt( tmp_U_amp4 + tmp_V_amp4 + 2tmp_U_amp2tmp_V_amp*2cos(2*(tmp_U_phi - tmp_V_phi)) )
1446	alpha2 = sqrt( tmpreal )
1447	if (alpha2 < 0) alpha2 = 0 !CYCLE
1448	alpha= sqrt( alpha2 )
1449
1450
1451	!major and minor axis of the ellipse
1452	qmax = sqrt( (tmp_U_amp2 + tmp_V_amp2 + alpha**2)/2 )
1453
1454	tmpreal = (tmp_U_amp2 + tmp_V_amp2 - alpha**2)/2
1455	if (tmpreal < 0) tmpreal = 0 !CYCLE
1456	!qmin = sqrt( (tmp_U_amp2 + tmp_V_amp2 - alpha**2)/2 ) ! but always positive.
1457	qmin = sqrt( tmpreal ) ! but always positive.
1458
1459	!eccentricity of ellipse
1460	tmpreal = (qmax + qmin)
1461	if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE
1462	!ecc = (qmax - qmin)/(qmax + qmin)
1463	ecc = (qmax - qmin)/(tmpreal)
1464
1465	! Angle of major and minor ellipse
1466	thetamax = atan2(( tmp_V_amp * cos((tmp_U_phi - tmp_V_phi) - delta) ) , ( tmp_U_amp * cos( delta) ) )
1467	thetamin = thetamax + rpi/2.
1468
1469
1470
1471	! Rotary current components: Pugh A3.10
1472	! Clockwise (c) and anticlockwise (ac) rotating rotate_wind_vectors
1473	! so Qc = clockwise = anticyclonic = negative
1474	! and Qac = anticlockwise = cyclonic = negative
1475
1476	tmpreal = tmp_U_amp2 + tmp_V_amp2 - (2tmp_U_amptmp_V_amp*sin( tmp_V_phi - tmp_U_phi))
1477	if (tmpreal < 0) tmpreal = 0 !CYCLE
1478	!Qc = 0.5sqrt( tmp_U_amp2 + tmp_V_amp2 - (2tmp_U_amptmp_V_ampsin( tmp_V_phi - tmp_U_phi)) )
1479	Qc = 0.5*sqrt( tmpreal )
1480
1481	tmpreal = tmp_U_amp2 + tmp_V_amp2 + (2tmp_U_amptmp_V_amp*sin( tmp_V_phi - tmp_U_phi))
1482	if (tmpreal < 0) tmpreal = 0 !CYCLE
1483	!Qac = 0.5sqrt( tmp_U_amp2 + tmp_V_amp2 + (2tmp_U_amptmp_V_ampsin( tmp_V_phi - tmp_U_phi)) )
1484	Qac = 0.5*sqrt( tmpreal )
1485
1486
1487	gc = atan2( ( ( tmp_U_ampsin( tmp_U_phi ) ) + (tmp_V_ampcos( tmp_V_phi) ) ) , &
1488	& ( (tmp_U_ampcos( tmp_U_phi )) - (tmp_V_ampsin( tmp_V_phi )) ) )
1489	gac = atan2( ( ( -tmp_U_ampsin( tmp_U_phi ) ) + (tmp_V_ampcos( tmp_V_phi) ) ) , &
1490	& ( (tmp_U_ampcos( tmp_U_phi )) + (tmp_V_ampsin( tmp_V_phi )) ) )
1491
1492	!Pugh A3.2
1493	Phi_Ua = -0.5*(gac - gc)
1494	dir_Ua = 0.5*(gac + gc) ! positive from x axis
1495
1496	tmpreal = qmax
1497	!if (tmpreal == 0) tmpreal = tmpreal + 0.0000001 !CYCLE
1498	if (tmpreal == 0) then
1499	polarity = 0
1500	else
1501	polarity = (Qac - Qc)/qmax
1502	endif
1503
1504
1505	tmp_u_amp_3d_mat(ji,jj,jk) = tmp_u_amp
1506	tmp_v_amp_3d_mat(ji,jj,jk) = tmp_v_amp
1507	tmp_u_phi_3d_mat(ji,jj,jk) = tmp_u_phi
1508	tmp_v_phi_3d_mat(ji,jj,jk) = tmp_v_phi
1509
1510
1511	a_u_3d_mat(ji,jj,jk) = a_u
1512	b_u_3d_mat(ji,jj,jk) = b_u
1513	a_v_3d_mat(ji,jj,jk) = a_v
1514	b_v_3d_mat(ji,jj,jk) = b_v
1515
1516	qmax_3d_mat(ji,jj,jk) = qmax
1517	qmin_3d_mat(ji,jj,jk) = qmin
1518
1519	ecc_3d_mat(ji,jj,jk) = ecc
1520	thetamax_3d_mat(ji,jj,jk) = thetamax
1521	thetamin_3d_mat(ji,jj,jk) = thetamin
1522
1523	Qc_3d_mat(ji,jj,jk) = Qc
1524	Qac_3d_mat(ji,jj,jk) = Qac
1525	gc_3d_mat(ji,jj,jk) = gc
1526	gac_3d_mat(ji,jj,jk) = gac
1527
1528	Phi_Ua_3d_mat(ji,jj,jk) = Phi_Ua
1529	dir_Ua_3d_mat(ji,jj,jk) = dir_Ua
1530	polarity_3d_mat(ji,jj,jk) = polarity
1531
1532	ENDIF
1533	END DO !ji
1534	END DO !jj
1535	END DO !jk
1536
1537
1538	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_amp_t_uv3d'
1539	IF( iom_use(TRIM(tmp_name)) ) THEN
1540	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1541	CALL iom_put( TRIM(tmp_name), tmp_u_amp_3d_mat(:,:,:))
1542	ENDIF
1543	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_amp_t_uv3d'
1544	IF( iom_use(TRIM(tmp_name)) ) THEN
1545	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1546	CALL iom_put( TRIM(tmp_name), tmp_v_amp_3d_mat(:,:,:))
1547	ENDIF
1548	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_u_phi_t_uv3d'
1549	IF( iom_use(TRIM(tmp_name)) ) THEN
1550	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1551	CALL iom_put( TRIM(tmp_name), tmp_u_phi_3d_mat(:,:,:))
1552	ENDIF
1553	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_v_phi_t_uv3d'
1554	IF( iom_use(TRIM(tmp_name)) ) THEN
1555	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1556	CALL iom_put( TRIM(tmp_name), tmp_v_phi_3d_mat(:,:,:))
1557	ENDIF
1558
1559
1560
1561	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_u_uv3d'
1562	IF( iom_use(TRIM(tmp_name)) ) THEN
1563	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1564	CALL iom_put( TRIM(tmp_name), a_u_3d_mat(:,:,:))
1565	ENDIF
1566	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_a_v_uv3d'
1567	IF( iom_use(TRIM(tmp_name)) ) THEN
1568	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1569	CALL iom_put( TRIM(tmp_name), a_v_3d_mat(:,:,:))
1570	ENDIF
1571	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_u_uv3d'
1572	IF( iom_use(TRIM(tmp_name)) ) THEN
1573	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1574	CALL iom_put( TRIM(tmp_name), b_u_3d_mat(:,:,:))
1575	ENDIF
1576	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_b_v_uv3d'
1577	IF( iom_use(TRIM(tmp_name)) ) THEN
1578	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1579	CALL iom_put( TRIM(tmp_name), b_v_3d_mat(:,:,:))
1580	ENDIF
1581
1582	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmax_uv3d'
1583	IF( iom_use(TRIM(tmp_name)) ) THEN
1584	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1585	CALL iom_put( TRIM(tmp_name), qmax_3d_mat(:,:,:))
1586	ENDIF
1587	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_qmin_uv3d'
1588	IF( iom_use(TRIM(tmp_name)) ) THEN
1589	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1590	CALL iom_put( TRIM(tmp_name), qmin_3d_mat(:,:,:))
1591	ENDIF
1592
1593	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_ecc_uv3d'
1594	IF( iom_use(TRIM(tmp_name)) ) THEN
1595	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1596	CALL iom_put( TRIM(tmp_name), ecc_3d_mat(:,:,:))
1597	ENDIF
1598	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamax_uv3d'
1599	IF( iom_use(TRIM(tmp_name)) ) THEN
1600	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1601	CALL iom_put( TRIM(tmp_name), thetamax_3d_mat(:,:,:))
1602	ENDIF
1603	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_thetamin_uv3d'
1604	IF( iom_use(TRIM(tmp_name)) ) THEN
1605	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1606	CALL iom_put( TRIM(tmp_name), thetamin_3d_mat(:,:,:))
1607	ENDIF
1608
1609	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qc_uv3d'
1610	IF( iom_use(TRIM(tmp_name)) ) THEN
1611	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1612	CALL iom_put( TRIM(tmp_name), Qc_3d_mat(:,:,:))
1613	ENDIF
1614	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Qac_uv3d'
1615	IF( iom_use(TRIM(tmp_name)) ) THEN
1616	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1617	CALL iom_put( TRIM(tmp_name), Qac_3d_mat(:,:,:))
1618	ENDIF
1619	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gc_uv3d'
1620	IF( iom_use(TRIM(tmp_name)) ) THEN
1621	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1622	CALL iom_put( TRIM(tmp_name), gc_3d_mat(:,:,:))
1623	ENDIF
1624	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_gac_uv3d'
1625	IF( iom_use(TRIM(tmp_name)) ) THEN
1626	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1627	CALL iom_put( TRIM(tmp_name), gac_3d_mat(:,:,:))
1628	ENDIF
1629
1630
1631	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_Phi_Ua_uv3d'
1632	IF( iom_use(TRIM(tmp_name)) ) THEN
1633	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1634	CALL iom_put( TRIM(tmp_name), Phi_Ua_3d_mat(:,:,:))
1635	ENDIF
1636	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_dir_Ua_uv3d'
1637	IF( iom_use(TRIM(tmp_name)) ) THEN
1638	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1639	CALL iom_put( TRIM(tmp_name), dir_Ua_3d_mat(:,:,:))
1640	ENDIF
1641	tmp_name=TRIM(Wave(ntide_all(jh))%cname_tide)//'_polarity_uv3d'
1642	IF( iom_use(TRIM(tmp_name)) ) THEN
1643	IF(lwp .AND. ln_diaharm_verbose) WRITE(numout,*) "diaharm_fast: iom_put: ",TRIM(tmp_name)
1644	CALL iom_put( TRIM(tmp_name), polarity_3d_mat(:,:,:))
1645	ENDIF
1646
1647	tmp_u_amp_3d_mat(:,:,:) = 0.
1648	tmp_v_amp_3d_mat(:,:,:) = 0.
1649	tmp_u_phi_3d_mat(:,:,:) = 0.
1650	tmp_v_phi_3d_mat(:,:,:) = 0.
1651
1652	a_u_3d_mat(:,:,:) = 0.
1653	b_u_3d_mat(:,:,:) = 0.
1654	a_v_3d_mat(:,:,:) = 0.
1655	b_v_3d_mat(:,:,:) = 0.
1656
1657	qmax_3d_mat(:,:,:) = 0.
1658	qmin_3d_mat(:,:,:) = 0.
1659
1660	ecc_3d_mat(:,:,:) = 0.
1661	thetamax_3d_mat(:,:,:) = 0.
1662	thetamin_3d_mat(:,:,:) = 0.
1663
1664	Qc_3d_mat(:,:,:) = 0.
1665	Qac_3d_mat(:,:,:) = 0.
1666	gc_3d_mat(:,:,:) = 0.
1667	gac_3d_mat(:,:,:) = 0.
1668
1669	Phi_Ua_3d_mat(:,:,:) = 0.
1670	dir_Ua_3d_mat(:,:,:) = 0.
1671	polarity_3d_mat(:,:,:) = 0.
1672
1673
1674	END DO !jh
1675
1676
1677
1678
1679
1680	IF(lwp) WRITE(numout,*) "diaharm_fast: Postprocess 3d velocity tidal parameters"
1681	ENDIF
1682
1683	CALL FLUSH(numout)
1684	ENDIF
1685
1686
1687	CALL FLUSH(numout)
1688
1689	! to output tidal parameters, u and v on t grid
1690	!
1691	! !== standard Cd ==!
1692	! DO jj = 2, jpjm1
1693	! DO ji = 2, jpim1
1694	! imk = k_mk(ji,jj) ! ocean bottom level at t-points
1695	! zut = un(ji,jj,imk) + un(ji-1,jj,imk) ! 2 x velocity at t-point
1696	! zvt = vn(ji,jj,imk) + vn(ji,jj-1,imk)
1697	! ! ! here pCd0 = maskboost drag
1698	! pCdU(ji,jj) = - pCd0(ji,jj) * SQRT( 0.25 * ( zutzut + zvtzvt ) + pke0 )
1699	! END DO
1700	! END DO
1701
1702
1703
1704	IF (ln_diaharm_postproc_vel) THEN
1705	IF (ln_ana_uvbar) THEN
1706
1707	DEALLOCATE(amp_u2d, amp_v2d, phi_u2d, phi_v2d )
1708
1709	DEALLOCATE(tmp_u_amp_2d_mat, tmp_v_amp_2d_mat, tmp_u_phi_2d_mat, tmp_v_phi_2d_mat )
1710
1711	DEALLOCATE(a_u_2d_mat, b_u_2d_mat, a_v_2d_mat, b_v_2d_mat )
1712	DEALLOCATE(qmax_2d_mat, qmin_2d_mat, ecc_2d_mat )
1713	DEALLOCATE(thetamax_2d_mat, thetamin_2d_mat, Qc_2d_mat, Qac_2d_mat)
1714	DEALLOCATE(gc_2d_mat, gac_2d_mat, Phi_Ua_2d_mat, dir_Ua_2d_mat)
1715	DEALLOCATE(polarity_2d_mat )
1716
1717	ENDIF
1718	IF(lwp) WRITE(numout,*) "diaharm_fast: Deallocated 2d velocity tidal parameters"
1719
1720	IF (ln_ana_uv3d) THEN
1721
1722	DEALLOCATE(amp_u3d, amp_v3d, phi_u3d, phi_v3d )
1723
1724	DEALLOCATE(tmp_u_amp_3d_mat, tmp_v_amp_3d_mat, tmp_u_phi_3d_mat, tmp_v_phi_3d_mat )
1725
1726	DEALLOCATE(a_u_3d_mat, b_u_3d_mat, a_v_3d_mat, b_v_3d_mat )
1727	DEALLOCATE(qmax_3d_mat, qmin_3d_mat, ecc_3d_mat )
1728	DEALLOCATE(thetamax_3d_mat, thetamin_3d_mat, Qc_3d_mat, Qac_3d_mat)
1729	DEALLOCATE(gc_3d_mat, gac_3d_mat, Phi_Ua_3d_mat, dir_Ua_3d_mat)
1730	DEALLOCATE(polarity_3d_mat )
1731
1732	ENDIF
1733	IF(lwp) WRITE(numout,*) "diaharm_fast: Deallocated 3d velocity tidal parameters"
1734
1735	ENDIF
1736
1737
1738	CALL FLUSH(numout)
1739	!
1740	END SUBROUTINE harm_ana_out
1741	!
1742	SUBROUTINE harm_rst_write(kt)
1743	!!----------------------------------------------------------------------
1744	!! * ROUTINE harm_ana_init *
1745	!!
1746	!! ** Purpose : To write out cummulated Tidal Harmomnic data to file for
1747	!! restarting
1748	!!
1749	!! ** Method : restart files will be dated by default
1750	!!
1751	!! ** input :
1752	!!
1753	!! ** Action : ...
1754	!!
1755	!! history :
1756	!! 0.0 ! 01-16 (Enda O'Dea) Original code
1757	!! ASSUMES dated file for rose , can change later to be more generic
1758	!!----------------------------------------------------------------------
1759	INTEGER, INTENT(in) :: kt ! ocean time-step
1760	!!
1761	INTEGER :: jh, j2d, j3d
1762	CHARACTER(LEN=20) :: clkt ! ocean time-step define as a character
1763	CHARACTER(LEN=50) :: clname ! ocean output restart file name
1764	CHARACTER(LEN=150) :: clpath ! full path to ocean output restart file
1765	CHARACTER(LEN=250) :: clfinal ! full name
1766
1767	!restart file
1768	DO j2d=1,nvar_2d
1769	CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
1770	DO jh=1,nb_ana
1771	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 ) )
1772	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 ) )
1773	ENDDO
1774	ENDDO
1775
1776	DO j3d=1,nvar_3d
1777	!JT CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
1778	CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName3d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
1779	DO jh=1,nb_ana
1780	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 ) )
1781	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 ) )
1782	ENDDO
1783	ENDDO
1784
1785	IF(lwp) THEN
1786	IF( kt > 999999999 ) THEN ; WRITE(clkt, * ) kt
1787	ELSE ; WRITE(clkt, '(i8.8)') kt
1788	ENDIF
1789	clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
1790	clpath = TRIM(cn_ocerst_outdir)
1791	IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
1792	IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for writing: ',TRIM(clpath)//clname
1793
1794	WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
1795	OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
1796	WRITE(66) cc
1797	WRITE(66) anau
1798	WRITE(66) anav
1799	WRITE(66) anaf
1800	WRITE(66) fjulday_startharm
1801	CLOSE(66)
1802	WRITE(numout,*) '----------------------------'
1803	WRITE(numout,*) ' harm_rst_write: DONE '
1804	WRITE(numout,*) cc
1805	WRITE(numout,*) anaf
1806	WRITE(numout,*) anau
1807	WRITE(numout,*) anav
1808	WRITE(numout,*) fjulday_startharm
1809	WRITE(numout,*) '----------------------------'
1810	ENDIF
1811
1812	END SUBROUTINE harm_rst_write
1813
1814	SUBROUTINE harm_rst_read
1815	!!----------------------------------------------------------------------
1816	!! * ROUTINE harm_ana_init *
1817	!!
1818	!! ** Purpose : To read in cummulated Tidal Harmomnic data to file for
1819	!! restarting
1820	!!
1821	!! ** Method :
1822	!!
1823	!! ** input :
1824	!!
1825	!! ** Action : ...
1826	!!
1827	!! history :
1828	!! 0.0 ! 01-16 (Enda O'Dea) Original code
1829	!! ASSUMES dated file for rose , can change later to be more generic
1830	!!----------------------------------------------------------------------
1831	CHARACTER(LEN=20) :: clkt ! ocean time-step define as a character
1832	CHARACTER(LEN=50) :: clname ! ocean output restart file name
1833	CHARACTER(LEN=150) :: clpath ! full path to ocean output restart file
1834	CHARACTER(LEN=250) :: clfinal ! full name
1835	INTEGER :: jh, j2d, j3d
1836
1837	IF( nit000 > 999999999 ) THEN ; WRITE(clkt, * ) nit000-1
1838	ELSE ; WRITE(clkt, '(i8.8)') nit000-1
1839	ENDIF
1840	clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
1841	clpath = TRIM(cn_ocerst_outdir)
1842	IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
1843
1844	IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for reading: ',TRIM(clpath)//clname
1845
1846	DO j2d=1,nvar_2d
1847	CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
1848	IF(lwp) WRITE(numout,*) "2D", j2d, m_posi_2d(j2d), m_varName2d( m_posi_2d(j2d) )
1849	DO jh=1,nb_ana
1850	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 ) )
1851	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 ) )
1852	ENDDO
1853	ENDDO
1854
1855	DO j3d=1,nvar_3d
1856	!JT CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
1857	CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName3d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
1858	IF(lwp) WRITE(numout,*) "3D", j3d, m_posi_3d(j3d), m_varName3d( m_posi_3d(j3d) )
1859
1860	DO jh=1,nb_ana
1861	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 ) )
1862	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 ) )
1863	ENDDO
1864	ENDDO
1865
1866	WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
1867	OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
1868	READ(66) cc
1869	READ(66) anau
1870	READ(66) anav
1871	READ(66) anaf
1872	READ(66) fjulday_startharm
1873	CLOSE(66)
1874
1875	IF(lwp) THEN
1876	WRITE(numout,*) '----------------------------'
1877	WRITE(numout,*) ' Checking anaf is correct'
1878	WRITE(numout,*) cc
1879	WRITE(numout,*) anaf
1880	WRITE(numout,*) fjulday_startharm
1881	WRITE(numout,*) '----------------------------'
1882	ENDIF
1883
1884	END SUBROUTINE harm_rst_read
1885
1886	!!======================================================================
1887
1888	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 @ 15608

Download in other formats: