Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

diaharm.F90

Last change on this file was 10685, checked in by jcastill, 5 years ago
Changes as Ash's files
File size: 29.8 KB

Line
1	MODULE diaharm
2	!!======================================================================
3	!! * MODULE diaharm *
4	!! Harmonic analysis of tidal constituents
5	!!======================================================================
6	!! History : 3.1 ! 2007 (O. Le Galloudec, J. Chanut) Original code
7	!!----------------------------------------------------------------------
8	#if defined key_diaharm
9	!!----------------------------------------------------------------------
10	!! 'key_diaharm'
11	!!
12	!! NB: 2017-12 : add 3D harmonic analysis of velocities
13	!! integration of Maria Luneva's development
14	!! 'key_3Ddiaharm'
15	!!----------------------------------------------------------------------
16	USE oce ! ocean dynamics and tracers variables
17	USE dom_oce ! ocean space and time domain
18	USE phycst
19	USE daymod
20	USE tide_mod
21	USE sbctide ! Tidal forcing or not
22	!
23	# if defined key_3Ddiaharm
24	USE zdf_oce
25	#endif
26	!
27	USE in_out_manager ! I/O units
28	USE iom ! I/0 library
29	USE ioipsl ! NetCDF IPSL library
30	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
31	USE timing ! preformance summary
32	USE wrk_nemo ! working arrays
33
34	IMPLICIT NONE
35	PRIVATE
36
37	LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .TRUE.
38
39	INTEGER, PARAMETER :: jpincomax = 2.*jpmax_harmo
40	INTEGER, PARAMETER :: jpdimsparse = jpincomax30024
41
42	! !! namelist variables
43	INTEGER :: nit000_han ! First time step used for harmonic analysis
44	INTEGER :: nitend_han ! Last time step used for harmonic analysis
45	INTEGER :: nstep_han ! Time step frequency for harmonic analysis
46	INTEGER :: nb_ana ! Number of harmonics to analyse
47
48
49	INTEGER , ALLOCATABLE, DIMENSION(:) :: name
50	REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, ut , vt , ft
51	# if defined key_3Ddiaharm
52	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:,:) :: ana_temp
53	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: out_eta , out_u, out_v , out_w , out_dzi
54	# else
55	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp
56	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta , out_u, out_v
57	# endif
58
59	INTEGER :: ninco, nsparse
60	INTEGER , DIMENSION(jpdimsparse) :: njsparse, nisparse
61	INTEGER , SAVE, DIMENSION(jpincomax) :: ipos1
62	REAL(wp), DIMENSION(jpdimsparse) :: valuesparse
63	REAL(wp), DIMENSION(jpincomax) :: ztmp4 , ztmp7
64	REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3 , zpilier
65	REAL(wp), SAVE, DIMENSION(jpincomax) :: zpivot
66
67	CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: tname ! Names of tidal constituents ('M2', 'K1',...)
68
69	PUBLIC dia_harm ! routine called by step.F90
70
71	!!----------------------------------------------------------------------
72	!! NEMO/OPA 3.5 , NEMO Consortium (2013)
73	!! $Id$
74	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
75	!!----------------------------------------------------------------------
76	CONTAINS
77
78	SUBROUTINE dia_harm_init
79	!!----------------------------------------------------------------------
80	!! * ROUTINE dia_harm_init *
81	!!
82	!! ** Purpose : Initialization of tidal harmonic analysis
83	!!
84	!! ** Method : Initialize frequency array and nodal factor for nit000_han
85	!!
86	!!--------------------------------------------------------------------
87	INTEGER :: jh, nhan, jl
88	INTEGER :: ios ! Local integer output status for namelist read
89
90	NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, tname
91	!!----------------------------------------------------------------------
92
93	IF(lwp) THEN
94	WRITE(numout,*)
95	WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization'
96	# if defined key_3Ddiaharm
97	WRITE(numout,*) ' - 3D harmonic analysis of currents actovated (key_3Ddiaharm)'
98	#endif
99	WRITE(numout,*) '~~~~~~~ '
100	ENDIF
101	!
102	IF( .NOT. ln_tide ) CALL ctl_stop( 'dia_harm_init : ln_tide must be true for harmonic analysis')
103	!
104	CALL tide_init_Wave
105	!
106	REWIND( numnam_ref ) ! Namelist nam_diaharm in reference namelist : Tidal harmonic analysis
107	READ ( numnam_ref, nam_diaharm, IOSTAT = ios, ERR = 901)
108	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm in reference namelist', lwp )
109
110	REWIND( numnam_cfg ) ! Namelist nam_diaharm in configuration namelist : Tidal harmonic analysis
111	READ ( numnam_cfg, nam_diaharm, IOSTAT = ios, ERR = 902 )
112	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm in configuration namelist', lwp )
113	IF(lwm) WRITE ( numond, nam_diaharm )
114	!
115	IF(lwp) THEN
116	WRITE(numout,*) 'First time step used for analysis: nit000_han= ', nit000_han
117	WRITE(numout,*) 'Last time step used for analysis: nitend_han= ', nitend_han
118	WRITE(numout,*) 'Time step frequency for harmonic analysis: nstep_han= ', nstep_han
119	ENDIF
120
121	! Basic checks on harmonic analysis time window:
122	! ----------------------------------------------
123	IF( nit000 > nit000_han ) CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000', &
124	& ' restart capability not implemented' )
125	IF( nitend < nitend_han ) CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend', &
126	& 'restart capability not implemented' )
127
128	IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 ) &
129	& CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' )
130
131	nb_ana = 0
132	DO jh=1,jpmax_harmo
133	DO jl=1,jpmax_harmo
134	IF(TRIM(tname(jh)) == Wave(jl)%cname_tide) THEN
135	nb_ana=nb_ana+1
136	ENDIF
137	END DO
138	END DO
139	!
140	IF(lwp) THEN
141	WRITE(numout,*) ' Namelist nam_diaharm'
142	WRITE(numout,*) ' nb_ana = ', nb_ana
143	CALL flush(numout)
144	ENDIF
145	!
146	IF (nb_ana > jpmax_harmo) THEN
147	IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nb_ana must be lower than jpmax_harmo, stop'
148	IF(lwp) WRITE(numout,*) ' jpmax_harmo= ', jpmax_harmo
149	nstop = nstop + 1
150	ENDIF
151
152	ALLOCATE(name (nb_ana))
153	DO jh=1,nb_ana
154	DO jl=1,jpmax_harmo
155	IF (TRIM(tname(jh)) .eq. Wave(jl)%cname_tide) THEN
156	name(jh) = jl
157	EXIT
158	END IF
159	END DO
160	END DO
161
162	! Initialize frequency array:
163	! ---------------------------
164	ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) )
165
166	CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana )
167
168	IF(lwp) WRITE(numout,*) 'Analysed frequency : ',nb_ana ,'Frequency '
169
170	DO jh = 1, nb_ana
171	IF(lwp) WRITE(numout,*) ' : ',tname(jh),' ',ana_freq(jh)
172	END DO
173
174	! Initialize temporary arrays:
175	! ----------------------------
176	# if defined key_3Ddiaharm
177	ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 5, jpk ) )
178	ana_temp(:,:,:,:,:) = 0._wp
179	# else
180	ALLOCATE( ana_temp( jpi, jpj, 2*nb_ana, 3 ) )
181	ana_temp(:,:,:,: ) = 0._wp
182	#endif
183
184	END SUBROUTINE dia_harm_init
185
186
187	SUBROUTINE dia_harm ( kt )
188	!!----------------------------------------------------------------------
189	!! * ROUTINE dia_harm *
190	!!
191	!! ** Purpose : Tidal harmonic analysis main routine
192	!!
193	!! ** Action : Sums ssh/u/v over time analysis [nit000_han,nitend_han]
194	!!
195	!!--------------------------------------------------------------------
196	INTEGER, INTENT( IN ) :: kt
197	!
198	INTEGER :: ji, jj, jh, jc, nhc
199	# if defined key_3Ddiaharm
200	INTEGER :: jk
201	# endif
202	REAL(wp) :: ztime, ztemp
203	!!--------------------------------------------------------------------
204	IF( nn_timing == 1 ) CALL timing_start('dia_harm')
205
206	IF( kt == nit000 ) CALL dia_harm_init
207
208	IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN
209
210	ztime = (kt-nit000+1) * rdt
211
212	nhc = 0
213	DO jh = 1, nb_ana
214	DO jc = 1, 2
215	nhc = nhc+1
216	ztemp =( MOD(jc,2) * ft(jh) COS(ana_freq(jh)ztime + vt(jh) + ut(jh)) &
217	& +(1.-MOD(jc,2))* ft(jh) SIN(ana_freq(jh)ztime + vt(jh) + ut(jh)))
218
219	! ssh, ub, vb are stored at the last level of 5d array
220	DO jj = 1,jpj
221	DO ji = 1,jpi
222	! Elevation and currents
223	# if defined key_3Ddiaharm
224	ana_temp(ji,jj,nhc,1,jpk) = ana_temp(ji,jj,nhc,1,jpk) + ztempsshn(ji,jj)ssmask (ji,jj)
225	ana_temp(ji,jj,nhc,2,jpk) = ana_temp(ji,jj,nhc,2,jpk) + ztempun_b(ji,jj)ssumask(ji,jj)
226	ana_temp(ji,jj,nhc,3,jpk) = ana_temp(ji,jj,nhc,3,jpk) + ztempvn_b(ji,jj)ssvmask(ji,jj)
227
228	ana_temp(ji,jj,nhc,5,jpk) = ana_temp(ji,jj,nhc,5,jpk) &
229	& + ztempbfrva(ji,jj)vn(ji,jj,mbkv(ji,jj))*ssvmask(ji,jj)
230	ana_temp(ji,jj,nhc,4,jpk) = ana_temp(ji,jj,nhc,4,jpk) &
231	& + ztempbfrua(ji,jj)un(ji,jj,mbku(ji,jj))*ssumask(ji,jj)
232	# else
233	ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztempsshn(ji,jj)ssmask (ji,jj)
234	ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztempun_b(ji,jj)ssumask(ji,jj)
235	ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztempvn_b(ji,jj)ssvmask(ji,jj)
236	# endif
237	END DO
238	END DO
239	!
240	# if defined key_3Ddiaharm
241	! 3d velocity and density:
242	DO jk=1,jpk-1
243	DO jj = 1,jpj
244	DO ji = 1,jpi
245	! density and velocity
246	ana_temp(ji,jj,nhc,1,jk) = ana_temp(ji,jj,nhc,1,jk) + ztemp*rhd(ji,jj,jk)
247	ana_temp(ji,jj,nhc,2,jk) = ana_temp(ji,jj,nhc,2,jk) + ztemp*(un(ji,jj,jk)-un_b(ji,jj)) &
248	& *umask(ji,jj,jk)
249	ana_temp(ji,jj,nhc,3,jk) = ana_temp(ji,jj,nhc,3,jk) + ztemp*(vn(ji,jj,jk)-vn_b(ji,jj)) &
250	& *vmask(ji,jj,jk)
251	ana_temp(ji,jj,nhc,4,jk) = ana_temp(ji,jj,nhc,4,jk) + ztemp*wn(ji,jj,jk)
252
253	ana_temp(ji,jj,nhc,5,jk) = ana_temp(ji,jj,nhc,5,jk) - 0.5gravztemp*(rhd(ji,jj,jk)+rhd(ji,jj,jk+1) )/max(rn2(ji,jj,jk),1.e-8_wp)
254	END DO
255	END DO
256	ENDDO
257	# endif
258
259	END DO
260	END DO
261	!
262	END IF
263
264	IF ( kt == nitend_han ) CALL dia_harm_end
265
266	IF( nn_timing == 1 ) CALL timing_stop('dia_harm')
267
268	END SUBROUTINE dia_harm
269
270
271	SUBROUTINE dia_harm_end
272	!!----------------------------------------------------------------------
273	!! * ROUTINE diaharm_end *
274	!!
275	!! ** Purpose : Compute the Real and Imaginary part of tidal constituents
276	!!
277	!! ** Action : Decompose the signal on the harmonic constituents
278	!!
279	!!--------------------------------------------------------------------
280	INTEGER :: ji, jj, jh, jc, jn, nhan, jl
281	# if defined key_3Ddiaharm
282	INTEGER :: jk
283	# endif
284	INTEGER :: ksp, kun, keq
285	REAL(wp) :: ztime, ztime_ini, ztime_end
286	REAL(wp) :: X1,X2
287	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ana_amp
288	!!--------------------------------------------------------------------
289	CALL wrk_alloc( jpi , jpj , jpmax_harmo , 2 , ana_amp )
290
291	IF(lwp) WRITE(numout,*)
292	IF(lwp) WRITE(numout,*) 'anharmo_end: kt=nitend_han: Perform harmonic analysis'
293	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
294
295	ztime_ini = nit000_han*rdt ! Initial time in seconds at the beginning of analysis
296	ztime_end = nitend_han*rdt ! Final time in seconds at the end of analysis
297	nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis
298
299	# if defined key_3Ddiaharm
300	ALLOCATE( out_eta(jpi,jpj,jpk,2*nb_ana), &
301	& out_u (jpi,jpj,jpk,2*nb_ana), &
302	& out_v (jpi,jpj,jpk,2*nb_ana), &
303	& out_w (jpi,jpj,jpk,2*nb_ana), &
304	& out_dzi(jpi,jpj,jpk,2*nb_ana) )
305	# else
306	ALLOCATE( out_eta(jpi,jpj,2*nb_ana), &
307	& out_u (jpi,jpj,2*nb_ana), &
308	& out_v (jpi,jpj,2*nb_ana) )
309	# endif
310
311	IF(lwp) WRITE(numout,*) 'ANA F OLD', ft
312	IF(lwp) WRITE(numout,*) 'ANA U OLD', ut
313	IF(lwp) WRITE(numout,*) 'ANA V OLD', vt
314
315
316	ninco = 2*nb_ana
317	ksp = 0
318	keq = 0
319	DO jn = 1, nhan
320	ztime=( (nhan-jn)ztime_ini + (jn-1)ztime_end )/FLOAT(nhan-1)
321	keq = keq + 1
322	kun = 0
323	DO jh = 1, nb_ana
324	DO jc = 1, 2
325	kun = kun + 1
326	ksp = ksp + 1
327	nisparse(ksp) = keq
328	njsparse(ksp) = kun
329	valuesparse(ksp) = ( MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &
330	& + (1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)) )
331	END DO
332	END DO
333	END DO
334
335	nsparse = ksp
336
337	! Density and Elevation:
338	# if defined key_3Ddiaharm
339	DO jk=1,jpk
340	# endif
341	DO jj = 1, jpj
342	DO ji = 1, jpi
343	! Fill input array
344	kun = 0
345	DO jh = 1, nb_ana
346	DO jc = 1, 2
347	kun = kun + 1
348	# if defined key_3Ddiaharm
349	ztmp4(kun)=ana_temp(ji,jj,kun,1,jk)
350	# else
351	ztmp4(kun)=ana_temp(ji,jj,kun,1)
352	# endif
353	END DO
354	END DO
355
356	CALL SUR_DETERMINE(jj)
357
358	! Fill output array
359	DO jh = 1, nb_ana
360	ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
361	ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
362	END DO
363	END DO
364	END DO
365
366
367	DO jj = 1, jpj
368	DO ji = 1, jpi
369	DO jh = 1, nb_ana
370	X1 = ana_amp(ji,jj,jh,1)
371	X2 =-ana_amp(ji,jj,jh,2)
372	# if defined key_3Ddiaharm
373	out_eta(ji,jj,jk,jh ) = X1 * tmask_i(ji,jj)
374	out_eta(ji,jj,jk,jh+nb_ana) = X2 * tmask_i(ji,jj)
375	# else
376	out_eta(ji,jj ,jh ) = X1 * tmask_i(ji,jj)
377	out_eta(ji,jj ,jh+nb_ana) = X2 * tmask_i(ji,jj)
378	# endif
379	END DO
380	END DO
381	END DO
382
383	! u-component of velocity
384	DO jj = 1, jpj
385	DO ji = 1, jpi
386	! Fill input array
387	kun=0
388	DO jh = 1,nb_ana
389	DO jc = 1,2
390	kun = kun + 1
391	# if defined key_3Ddiaharm
392	ztmp4(kun)=ana_temp(ji,jj,kun,2,jk)
393	# else
394	ztmp4(kun)=ana_temp(ji,jj,kun,2)
395	# endif
396	END DO
397	END DO
398
399	CALL SUR_DETERMINE(jj+1)
400
401	! Fill output array
402	DO jh = 1, nb_ana
403	ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1)
404	ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2)
405	END DO
406
407	END DO
408	END DO
409
410	DO jj = 1, jpj
411	DO ji = 1, jpi
412	DO jh = 1, nb_ana
413	X1= ana_amp(ji,jj,jh,1)
414	X2=-ana_amp(ji,jj,jh,2)
415	# if defined key_3Ddiaharm
416	out_u(ji,jj,jk, jh) = X1 * ssumask(ji,jj)
417	out_u(ji,jj,jk,nb_ana+jh) = X2 * ssumask(ji,jj)
418	# else
419	out_u(ji,jj, jh) = X1 * ssumask(ji,jj)
420	out_u(ji,jj, nb_ana+jh) = X2 * ssumask(ji,jj)
421	# endif
422	ENDDO
423	ENDDO
424	ENDDO
425
426	! v- velocity
427	DO jj = 1, jpj
428	DO ji = 1, jpi
429	! Fill input array
430	kun=0
431	DO jh = 1,nb_ana
432	DO jc = 1,2
433	kun = kun + 1
434	# if defined key_3Ddiaharm
435	ztmp4(kun)=ana_temp(ji,jj,kun,3,jk)
436	# else
437	ztmp4(kun)=ana_temp(ji,jj,kun,3)
438	# endif
439	END DO
440	END DO
441
442	CALL SUR_DETERMINE(jj+1)
443
444	! Fill output array
445	DO jh = 1, nb_ana
446	ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
447	ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
448	END DO
449
450	END DO
451	END DO
452
453	DO jj = 1, jpj
454	DO ji = 1, jpi
455	DO jh = 1, nb_ana
456	X1=ana_amp(ji,jj,jh,1)
457	X2=-ana_amp(ji,jj,jh,2)
458	# if defined key_3Ddiaharm
459	out_v(ji,jj,jk, jh)=X1 * ssvmask(ji,jj)
460	out_v(ji,jj,jk,nb_ana+jh)=X2 * ssvmask(ji,jj)
461	# else
462	out_v(ji,jj, jh)=X1 * ssvmask(ji,jj)
463	out_v(ji,jj, nb_ana+jh)=X2 * ssvmask(ji,jj)
464	# endif
465	END DO
466	END DO
467	END DO
468
469	# if defined key_3Ddiaharm
470	! w- velocity
471	DO jj = 1, jpj
472	DO ji = 1, jpi
473	! Fill input array
474	kun=0
475	DO jh = 1,nb_ana
476	DO jc = 1,2
477	kun = kun + 1
478	ztmp4(kun)=ana_temp(ji,jj,kun,4,jk)
479	END DO
480	END DO
481
482	CALL SUR_DETERMINE(jj+1)
483
484	! Fill output array
485	DO jh = 1, nb_ana
486	ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
487	ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
488	END DO
489
490	END DO
491	END DO
492
493	DO jj = 1, jpj
494	DO ji = 1, jpi
495	DO jh = 1, nb_ana
496	X1=ana_amp(ji,jj,jh,1)
497	X2=-ana_amp(ji,jj,jh,2)
498	out_w(ji,jj,jk, jh)=X1 * tmask_i(ji,jj)
499	out_w(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj)
500	END DO
501	END DO
502	END DO
503
504	! dzi- isopycnal displacements
505	DO jj = 1, jpj
506	DO ji = 1, jpi
507	! Fill input array
508	kun=0
509	DO jh = 1,nb_ana
510	DO jc = 1,2
511	kun = kun + 1
512	ztmp4(kun)=ana_temp(ji,jj,kun,5,jk)
513	END DO
514	END DO
515
516	CALL SUR_DETERMINE(jj+1)
517
518	! Fill output array
519	DO jh = 1, nb_ana
520	ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
521	ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
522	END DO
523
524	END DO
525	END DO
526
527	DO jj = 1, jpj
528	DO ji = 1, jpi
529	DO jh = 1, nb_ana
530	X1=ana_amp(ji,jj,jh,1)
531	X2=-ana_amp(ji,jj,jh,2)
532	out_dzi(ji,jj,jk, jh)=X1 * tmask_i(ji,jj)
533	out_dzi(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj)
534	END DO
535	END DO
536	END DO
537
538	ENDDO ! jk
539	# endif
540
541	CALL dia_wri_harm ! Write results in files
542	CALL wrk_dealloc( jpi , jpj , jpmax_harmo , 2 , ana_amp )
543	!
544	END SUBROUTINE dia_harm_end
545
546
547	SUBROUTINE dia_wri_harm
548	!!--------------------------------------------------------------------
549	!! * ROUTINE dia_wri_harm *
550	!!
551	!! ** Purpose : Write tidal harmonic analysis results in a netcdf file
552	!!--------------------------------------------------------------------
553	CHARACTER(LEN=lc) :: cltext
554	CHARACTER(LEN=lc) :: &
555	cdfile_name_T , & ! name of the file created (T-points)
556	cdfile_name_U , & ! name of the file created (U-points)
557	cdfile_name_V ! name of the file created (V-points)
558	INTEGER :: jh
559
560	# if defined key_3Ddiaharm
561	CHARACTER(LEN=lc) :: cdfile_name_W ! name of the file created (W-points)
562	INTEGER :: jk
563	REAL(WP), ALLOCATABLE, DIMENSION (:,:,:) :: z3real, z3im
564	REAL(WP), ALLOCATABLE, DIMENSION (:,:) :: z2real, z2im
565	# endif
566	!!----------------------------------------------------------------------
567
568	#if defined key_dimgout
569	cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc'
570	cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc'
571	cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc'
572	# if defined key_3Ddiaharm
573	cdfile_name_W = TRIM(cexper)//'_Tidal_harmonics_gridW.dimgproc'
574	# endif
575	#endif
576
577	IF(lwp) WRITE(numout,*) ' '
578	IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
579	#if defined key_dimgout
580	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ Output files: ', TRIM(cdfile_name_T)
581	IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_U)
582	IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_V)
583	# if defined key_3Ddiaharm
584	IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_W)
585	# endif
586	#endif
587	IF(lwp) WRITE(numout,*) ' '
588
589	# if defined key_3Ddiaharm
590	ALLOCATE( z3real(jpi,jpj,jpk),z3im(jpi,jpj,jpk),z2real(jpi,jpj),z2im(jpi,jpj))
591	# endif
592
593	! A) density and elevation
594	!/////////////
595	!
596	#if defined key_dimgout
597	cltext='density amplitude and phase; elevation is level=jpk '
598	CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2')
599	#else
600	# if defined key_3Ddiaharm
601	z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp
602	# endif
603	DO jh = 1, nb_ana
604	# if defined key_3Ddiaharm
605	DO jk=1,jpkm1
606	z3real(:,:,jk)=out_eta(:,:,jk,jh)
607	z3im (:,:,jk)=out_eta(:,:,jk,jh+nb_ana)
608	ENDDO
609	z2real(:,:)=out_eta(:,:,jpk,jh); z2im(:,:)=out_eta(:,:,jpk,jh+nb_ana)
610	CALL iom_put( TRIM(tname(jh))//'x_ro', z3real(:,:,:) )
611	CALL iom_put( TRIM(tname(jh))//'y_ro', z3im (:,:,:) )
612	CALL iom_put( TRIM(tname(jh))//'x' , z2real(:,: ) )
613	CALL iom_put( TRIM(tname(jh))//'y' , z2im (:,: ) )
614	# else
615	WRITE(numout,*) "OUTPUT ORI: ", TRIM(tname(jh))//'x', ' & ', TRIM(tname(jh))//'y', MAXVAL(out_eta(:,:,jh))
616	CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
617	CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) )
618	# endif
619	END DO
620	#endif
621
622	! B) u
623	!/////////
624	!
625	#if defined key_dimgout
626	cltext='3d u amplitude and phase; ubar is the last level'
627	CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2')
628	#else
629	# if defined key_3Ddiaharm
630	z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp
631	# endif
632	DO jh = 1, nb_ana
633	# if defined key_3Ddiaharm
634	DO jk=1,jpkm1
635	z3real(:,:,jk)=out_u(:,:,jk,jh)
636	z3im (:,:,jk)=out_u(:,:,jk,jh+nb_ana)
637	ENDDO
638	z2real(:,:)=out_u(:,:,jpk,jh); z2im(:,:)=out_u(:,:,jpk,jh+nb_ana)
639	CALL iom_put( TRIM(tname(jh))//'x_u3d', z3real(:,:,:) )
640	CALL iom_put( TRIM(tname(jh))//'y_u3d', z3im (:,:,:) )
641	CALL iom_put( TRIM(tname(jh))//'x_u2d', z2real(:,:) )
642	CALL iom_put( TRIM(tname(jh))//'y_u2d', z2im (:,:) )
643	z2real(:,:)=out_w(:,:,jpk,jh); z2im(:,:)=out_w(:,:,jpk,jh+nb_ana)
644	CALL iom_put( TRIM(tname(jh))//'x_tabx', z2real(:,:) )
645	CALL iom_put( TRIM(tname(jh))//'y_tabx', z2im (:,:) )
646	# else
647	CALL iom_put( TRIM(tname(jh))//'x_u2d', out_u(:,:,jh) )
648	CALL iom_put( TRIM(tname(jh))//'y_u2d', out_u(:,:,nb_ana+jh) )
649	# endif
650	END DO
651	#endif
652
653	! C) v
654	!/////////
655	!
656	#if defined key_dimgout
657	cltext='3d v amplitude and phase; vbar is the last level'
658	CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2')
659	#else
660	# if defined key_3Ddiaharm
661	z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp
662	# endif
663	DO jh = 1, nb_ana
664	# if defined key_3Ddiaharm
665	DO jk=1,jpkm1
666	z3real(:,:,jk)=out_v(:,:,jk,jh)
667	z3im (:,:,jk)=out_v(:,:,jk,jh+nb_ana)
668	ENDDO
669	z2real(:,:)=out_v(:,:,jpk,jh); z2im(:,:)=out_v(:,:,jpk,jh+nb_ana)
670	CALL iom_put( TRIM(tname(jh))//'x_v3d', z3real(:,:,:) )
671	CALL iom_put( TRIM(tname(jh))//'y_v3d', z3im (:,:,:) )
672	CALL iom_put( TRIM(tname(jh))//'x_v2d' , z2real(:,:) )
673	CALL iom_put( TRIM(tname(jh))//'y_v2d' , z2im (:,:) )
674	z2real(:,:)=out_dzi(:,:,jpk,jh); z2im(:,:)=out_dzi(:,:,jpk,jh+nb_ana)
675	CALL iom_put( TRIM(tname(jh))//'x_taby', z2real(:,:) )
676	CALL iom_put( TRIM(tname(jh))//'y_taby', z2im (:,:) )
677	# else
678	CALL iom_put( TRIM(tname(jh))//'x_v2d', out_v(:,:,jh ) )
679	CALL iom_put( TRIM(tname(jh))//'y_v2d', out_v(:,:,jh+nb_ana) )
680	# endif
681	END DO
682
683	#endif
684	! D) w
685	# if defined key_3Ddiaharm
686	# if defined key_dimgout
687	cltext='3d w amplitude and phase; vort_baro is the last level'
688	CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2')
689	# else
690	DO jh = 1, nb_ana
691	DO jk=1,jpkm1
692	z3real(:,:,jk)=out_w(:,:,jk,jh)
693	z3im(:,:,jk)=out_w(:,:,jk,jh+nb_ana)
694	ENDDO
695	CALL iom_put( TRIM(tname(jh))//'x_w3d', z3real(:,:,:) )
696	CALL iom_put( TRIM(tname(jh))//'y_w3d', z3im(:,:,:) )
697	END DO
698	# endif
699
700	! E) dzi + tau_bot
701	# if defined key_dimgout
702	cltext='dzi=g*ro/N2 amplitude and phase'
703	CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2')
704	# else
705	DO jh = 1, nb_ana
706	DO jk=1,jpkm1
707	z3real(:,:,jk)=out_dzi(:,:,jk,jh)
708	z3im(:,:,jk)=out_dzi(:,:,jk,jh+nb_ana)
709	ENDDO
710	CALL iom_put( TRIM(tname(jh))//'x_dzi', z3real(:,:,:) )
711	CALL iom_put( TRIM(tname(jh))//'y_dzi', z3im(:,:,:) )
712	END DO
713	# endif
714	# endif
715
716	!
717	# if defined key_3Ddiaharm
718	DEALLOCATE(z3real, z3im, z2real,z2im)
719	# endif
720
721	END SUBROUTINE dia_wri_harm
722
723
724	SUBROUTINE SUR_DETERMINE(init)
725	!!---------------------------------------------------------------------------------
726	!! * ROUTINE SUR_DETERMINE *
727	!!
728	!!
729	!!
730	!!---------------------------------------------------------------------------------
731	INTEGER, INTENT(in) :: init
732	!
733	INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
734	REAL(wp) :: zval1, zval2, zx1
735	REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2
736	INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot
737	!---------------------------------------------------------------------------------
738	CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 )
739	CALL wrk_alloc( jpincomax , ipos2 , ipivot )
740
741	IF( init == 1 ) THEN
742	IF( nsparse > jpdimsparse ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
743	IF( ninco > jpincomax ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
744	!
745	ztmp3(:,:) = 0._wp
746	!
747	DO jk1_sd = 1, nsparse
748	DO jk2_sd = 1, nsparse
749	nisparse(jk2_sd) = nisparse(jk2_sd)
750	njsparse(jk2_sd) = njsparse(jk2_sd)
751	IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
752	ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) &
753	& + valuesparse(jk1_sd)*valuesparse(jk2_sd)
754	ENDIF
755	END DO
756	END DO
757	!
758	DO jj_sd = 1 ,ninco
759	ipos1(jj_sd) = jj_sd
760	ipos2(jj_sd) = jj_sd
761	END DO
762	!
763	DO ji_sd = 1 , ninco
764	!
765	!find greatest non-zero pivot:
766	zval1 = ABS(ztmp3(ji_sd,ji_sd))
767	!
768	ipivot(ji_sd) = ji_sd
769	DO jj_sd = ji_sd, ninco
770	zval2 = ABS(ztmp3(ji_sd,jj_sd))
771	IF( zval2 >= zval1 )THEN
772	ipivot(ji_sd) = jj_sd
773	zval1 = zval2
774	ENDIF
775	END DO
776	!
777	DO ji1_sd = 1, ninco
778	zcol1(ji1_sd) = ztmp3(ji1_sd,ji_sd)
779	zcol2(ji1_sd) = ztmp3(ji1_sd,ipivot(ji_sd))
780	ztmp3(ji1_sd,ji_sd) = zcol2(ji1_sd)
781	ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)
782	END DO
783	!
784	ipos2(ji_sd) = ipos1(ipivot(ji_sd))
785	ipos2(ipivot(ji_sd)) = ipos1(ji_sd)
786	ipos1(ji_sd) = ipos2(ji_sd)
787	ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))
788	zpivot(ji_sd) = ztmp3(ji_sd,ji_sd)
789	DO jj_sd = 1, ninco
790	ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)
791	END DO
792	!
793	DO ji2_sd = ji_sd+1, ninco
794	zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
795	DO jj_sd=1,ninco
796	ztmp3(ji2_sd,jj_sd)= ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
797	END DO
798	END DO
799	!
800	END DO
801	!
802	ENDIF ! End init==1
803
804	DO ji_sd = 1, ninco
805	ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)
806	DO ji2_sd = ji_sd+1, ninco
807	ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
808	END DO
809	END DO
810
811	!system solving:
812	ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco)
813	ji_sd = ninco
814	DO ji_sd = ninco-1, 1, -1
815	zx1 = 0._wp
816	DO jj_sd = ji_sd+1, ninco
817	zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
818	END DO
819	ztmpx(ji_sd) = ztmp4(ji_sd)-zx1
820	END DO
821
822	DO jj_sd =1, ninco
823	ztmp7(ipos1(jj_sd))=ztmpx(jj_sd)
824	END DO
825
826	CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 )
827	CALL wrk_dealloc( jpincomax , ipos2 , ipivot )
828	!
829	END SUBROUTINE SUR_DETERMINE
830
831	#else
832	!!----------------------------------------------------------------------
833	!! Default case : Empty module
834	!!----------------------------------------------------------------------
835	LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .FALSE.
836	CONTAINS
837	SUBROUTINE dia_harm ( kt ) ! Empty routine
838	INTEGER, INTENT( IN ) :: kt
839	WRITE(,) 'dia_harm: you should not have seen this print'
840	END SUBROUTINE dia_harm
841	#endif
842
843	!!======================================================================
844	END MODULE diaharm

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

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

Context Navigation

source: branches/UKMO/r8395_India_uncoupled/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90

Download in other formats: