Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

diaharm.F90 @ 12453

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

Download in other formats: