Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

diaharm.F90 @ 15422

Last change on this file since 15422 was 15422, checked in by jcastill, 3 years ago
Changes tested so that they can merged with the CO9 Met Office branch - jpmax_harmo should be 34 with FES14 tides, but the last components are not used anyway
File size: 28.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	!! 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 = jpincomax36624*2 ! 30min for a 1yr-long run
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	REAL(wp) :: z1_tmp3
55	INTEGER , DIMENSION(jpdimsparse) :: njsparse, nisparse
56	INTEGER , SAVE, DIMENSION(jpincomax) :: ipos1
57	REAL(wp), DIMENSION(jpdimsparse) :: valuesparse
58	REAL(wp), DIMENSION(jpincomax) :: ztmp4 , ztmp7, z1_pivot
59	REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3 , zpilier
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 = 2, jpjm1
218	DO ji = 2, jpim1
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	END IF
258	!
259	IF( kt == nitend_han ) CALL dia_harm_end
260	!
261	IF( ln_timing ) CALL timing_stop('dia_harm')
262	!
263	END SUBROUTINE dia_harm
264
265
266	SUBROUTINE dia_harm_end
267	!!----------------------------------------------------------------------
268	!! * ROUTINE diaharm_end *
269	!!
270	!! ** Purpose : Compute the Real and Imaginary part of tidal constituents
271	!!
272	!! ** Action : Decompose the signal on the harmonic constituents
273	!!
274	!!--------------------------------------------------------------------
275	INTEGER :: ji, jj, jh, jc, jn, nhan
276	# if defined key_3Ddiaharm
277	INTEGER :: jk
278	# endif
279	INTEGER :: ksp, kun, keq
280	REAL(wp) :: ztime, ztime_ini, ztime_end, z1_han
281	!!--------------------------------------------------------------------
282	!
283	IF(lwp) WRITE(numout,*)
284	IF(lwp) WRITE(numout,*) 'dia_harm_end: kt=nitend_han: Perform harmonic analysis'
285	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
286
287	ztime_ini = nit000_han*rdt ! Initial time in seconds at the beginning of analysis
288	ztime_end = nitend_han*rdt ! Final time in seconds at the end of analysis
289	nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis
290	z1_han = 1._wp / REAL(nhan-1)
291
292	# if defined key_3Ddiaharm
293	ALLOCATE( out_eta(jpi,jpj,jpk,2*nb_ana), &
294	& out_u (jpi,jpj,jpk,2*nb_ana), &
295	& out_v (jpi,jpj,jpk,2*nb_ana), &
296	& out_w (jpi,jpj,jpk,2*nb_ana), &
297	& out_dzi(jpi,jpj,jpk,2*nb_ana) )
298	# else
299	ALLOCATE( out_eta(jpi,jpj,2*nb_ana), &
300	& out_u (jpi,jpj,2*nb_ana), &
301	& out_v (jpi,jpj,2*nb_ana) )
302	# endif
303
304	IF(lwp) WRITE(numout,*) 'ANA F OLD', ft
305	IF(lwp) WRITE(numout,*) 'ANA U OLD', ut
306	IF(lwp) WRITE(numout,*) 'ANA V OLD', vt
307
308	ninco = 2*nb_ana
309
310	ksp = 0
311	keq = 0
312	DO jn = 1, nhan
313	ztime=( (nhan-jn)ztime_ini + (jn-1)ztime_end ) * z1_han
314	keq = keq + 1
315	kun = 0
316	DO jh = 1, nb_ana
317	DO jc = 1, 2
318	kun = kun + 1
319	ksp = ksp + 1
320	nisparse(ksp) = keq
321	njsparse(ksp) = kun
322	valuesparse(ksp) = ( MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &
323	& + (1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)) )
324	END DO
325	END DO
326	END DO
327
328	nsparse = ksp
329
330	IF( nsparse > jpdimsparse ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
331	IF( ninco > jpincomax ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
332
333	CALL SUR_DETERMINE_INIT
334
335	! Density and Elevation:
336	# if defined key_3Ddiaharm
337	DO jk=1,jpk
338	# endif
339	DO jj = 2, jpjm1
340	DO ji = 2, jpim1
341
342	! Fill input array
343	# if defined key_3Ddiaharm
344	ztmp4(1:nb_ana2) = ana_temp(ji,jj,1:nb_ana2,1,jk)
345	# else
346	ztmp4(1:nb_ana2) = ana_temp(ji,jj,1:nb_ana2,1)
347	# endif
348	CALL SUR_DETERMINE
349
350	! Fill output array
351	DO jh = 1, nb_ana
352	# if defined key_3Ddiaharm
353	out_eta(ji,jj,jk,jh ) = ztmp7((jh-1)2+1) ssmask(ji,jj)
354	out_eta(ji,jj,jk,jh+nb_ana) = -ztmp7((jh-1)2+2) ssmask(ji,jj)
355	# else
356	out_eta(ji,jj, jh ) = ztmp7((jh-1)2+1) ssmask(ji,jj)
357	out_eta(ji,jj, jh+nb_ana) = -ztmp7((jh-1)2+2) ssmask(ji,jj)
358	# endif
359	END DO
360	END DO
361	END DO
362
363	! ubar:
364	DO jj = 2, jpjm1
365	DO ji = 2, jpim1
366
367	! Fill input array
368	# if defined key_3Ddiaharm
369	ztmp4(1:nb_ana2) = ana_temp(ji,jj,1:nb_ana2,2,jk)
370	# else
371	ztmp4(1:nb_ana2) = ana_temp(ji,jj,1:nb_ana2,2)
372	# endif
373	CALL SUR_DETERMINE
374
375	! Fill output array
376	DO jh = 1, nb_ana
377	# if defined key_3Ddiaharm
378	out_u(ji,jj,jk, jh) = ztmp7((jh-1)2+1) ssumask(ji,jj)
379	out_u(ji,jj,jk,nb_ana+jh) = -ztmp7((jh-1)2+2) ssumask(ji,jj)
380	# else
381	out_u(ji,jj, jh) = ztmp7((jh-1)2+1) ssumask(ji,jj)
382	out_u(ji,jj, nb_ana+jh) = -ztmp7((jh-1)2+2) ssumask(ji,jj)
383	# endif
384
385	END DO
386
387	END DO
388	END DO
389
390	! vbar:
391	DO jj = 2, jpjm1
392	DO ji = 2, jpim1
393
394	! Fill input array
395	# if defined key_3Ddiaharm
396	ztmp4(1:nb_ana2) = ana_temp(ji,jj,1:nb_ana2,3,jk)
397	# else
398	ztmp4(1:nb_ana2) = ana_temp(ji,jj,1:nb_ana2,3)
399	# endif
400	CALL SUR_DETERMINE
401
402	! Fill output array
403	DO jh = 1, nb_ana
404	# if defined key_3Ddiaharm
405	out_v(ji,jj,jk, jh) = ztmp7((jh-1)2+1) ssvmask(ji,jj)
406	out_v(ji,jj,jk,nb_ana+jh) = -ztmp7((jh-1)2+2) ssvmask(ji,jj)
407	# else
408	out_v(ji,jj, jh) = ztmp7((jh-1)2+1) ssvmask(ji,jj)
409	out_v(ji,jj, nb_ana+jh) = -ztmp7((jh-1)2+2) ssvmask(ji,jj)
410	# endif
411	END DO
412
413	END DO
414	END DO
415
416	# if defined key_3Ddiaharm
417	! w- velocity
418	DO jj = 1, jpj
419	DO ji = 1, jpi
420	! Fill input array
421	kun=0
422	DO jh = 1,nb_ana
423	DO jc = 1,2
424	kun = kun + 1
425	ztmp4(kun)=ana_temp(ji,jj,kun,4,jk)
426	END DO
427	END DO
428
429	CALL SUR_DETERMINE(jj+1)
430
431	! Fill output array
432	DO jh = 1, nb_ana
433	ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
434	ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
435	END DO
436
437	END DO
438	END DO
439
440	DO jj = 1, jpj
441	DO ji = 1, jpi
442	DO jh = 1, nb_ana
443	X1=ana_amp(ji,jj,jh,1)
444	X2=-ana_amp(ji,jj,jh,2)
445	out_w(ji,jj,jk, jh)=X1 * tmask_i(ji,jj)
446	out_w(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj)
447	END DO
448	END DO
449	END DO
450
451	! dzi- isopycnal displacements
452	DO jj = 1, jpj
453	DO ji = 1, jpi
454	! Fill input array
455	kun=0
456	DO jh = 1,nb_ana
457	DO jc = 1,2
458	kun = kun + 1
459	ztmp4(kun)=ana_temp(ji,jj,kun,5,jk)
460	END DO
461	END DO
462
463	CALL SUR_DETERMINE(jj+1)
464
465	! Fill output array
466	DO jh = 1, nb_ana
467	ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
468	ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
469	END DO
470
471	END DO
472	END DO
473
474	DO jj = 1, jpj
475	DO ji = 1, jpi
476	DO jh = 1, nb_ana
477	X1=ana_amp(ji,jj,jh,1)
478	X2=-ana_amp(ji,jj,jh,2)
479	out_dzi(ji,jj,jk, jh)=X1 * tmask_i(ji,jj)
480	out_dzi(ji,jj,jk,nb_ana+jh)=X2 * tmask_i(ji,jj)
481	END DO
482	END DO
483	END DO
484
485	ENDDO ! jk
486	# endif
487	!
488	! clem: we could avoid this call if all the loops were from 1:jpi and 1:jpj
489	! but I think this is the most efficient
490	CALL lbc_lnk_multi( 'dia_harm_end', out_eta, 'T', 1., out_u, 'U', -1. , out_v, 'V', -1. )
491	!
492	CALL dia_wri_harm ! Write results in files
493	!
494	DEALLOCATE( out_eta, out_u, out_v )
495	!
496	END SUBROUTINE dia_harm_end
497
498
499	SUBROUTINE dia_wri_harm
500	!!--------------------------------------------------------------------
501	!! * ROUTINE dia_wri_harm *
502	!!
503	!! ** Purpose : Write tidal harmonic analysis results in a netcdf file
504	!!--------------------------------------------------------------------
505	INTEGER :: jh
506
507	# if defined key_3Ddiaharm
508	CHARACTER(LEN=lc) :: cdfile_name_W ! name of the file created (W-points)
509	INTEGER :: jk
510	REAL(WP), ALLOCATABLE, DIMENSION (:,:,:) :: z3real, z3im
511	REAL(WP), ALLOCATABLE, DIMENSION (:,:) :: z2real, z2im
512	# endif
513	!!----------------------------------------------------------------------
514
515	#if defined key_dimgout
516	cdfile_name_T = TRIM(cexper)//'_Tidal_harmonics_gridT.dimgproc'
517	cdfile_name_U = TRIM(cexper)//'_Tidal_harmonics_gridU.dimgproc'
518	cdfile_name_V = TRIM(cexper)//'_Tidal_harmonics_gridV.dimgproc'
519	# if defined key_3Ddiaharm
520	cdfile_name_W = TRIM(cexper)//'_Tidal_harmonics_gridW.dimgproc'
521	# endif
522	#endif
523
524	IF(lwp) WRITE(numout,*) ' '
525	IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
526	#if defined key_dimgout
527	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ Output files: ', TRIM(cdfile_name_T)
528	IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_U)
529	IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_V)
530	# if defined key_3Ddiaharm
531	IF(lwp) WRITE(numout,*) ' ', TRIM(cdfile_name_W)
532	# endif
533	#endif
534	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
535
536	# if defined key_3Ddiaharm
537	ALLOCATE(z3real(jpi,jpj,jpk),z3im(jpi,jpj,jpk),z2real(jpi,jpj),z2im(jpi,jpj))
538	# endif
539
540	! A) density and elevation
541	!/////////////
542	#if defined key_dimgout
543	cltext='density amplitude and phase; elevation is level=jpk '
544	CALL dia_wri_dimg(TRIM(cdfile_name_T), TRIM(cltext), out_eta, 2*nb_ana, '2')
545	#else
546	# if defined key_3Ddiaharm
547	z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp
548	# endif
549	DO jh = 1, nb_ana
550	# if defined key_3Ddiaharm
551	DO jk=1,jpkm1
552	z3real(:,:,jk)=out_eta(:,:,jk,jh)
553	z3im (:,:,jk)=out_eta(:,:,jk,jh+nb_ana)
554	ENDDO
555	z2real(:,:)=out_eta(:,:,jpk,jh); z2im(:,:)=out_eta(:,:,jpk,jh+nb_ana)
556	CALL iom_put( TRIM(tname(jh))//'x_ro', z3real(:,:,:) )
557	CALL iom_put( TRIM(tname(jh))//'y_ro', z3im (:,:,:) )
558	CALL iom_put( TRIM(tname(jh))//'x' , z2real(:,: ) )
559	CALL iom_put( TRIM(tname(jh))//'y' , z2im (:,: ) )
560	# else
561	WRITE(numout,*) "OUTPUT ORI: ", TRIM(tname(jh))//'x', ' & ', TRIM(tname(jh))//'y', MAXVAL(out_eta(:,:,jh))
562	CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
563	CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,jh+nb_ana) )
564	# endif
565	END DO
566	#endif
567
568	! B) u
569	!/////////
570	#if defined key_dimgout
571	cltext='3d u amplitude and phase; ubar is the last level'
572	CALL dia_wri_dimg(TRIM(cdfile_name_U), TRIM(cltext), out_u, 2*nb_ana, '2')
573	#else
574	# if defined key_3Ddiaharm
575	z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp
576	# endif
577	DO jh = 1, nb_ana
578	# if defined key_3Ddiaharm
579	DO jk=1,jpkm1
580	z3real(:,:,jk)=out_u(:,:,jk,jh)
581	z3im (:,:,jk)=out_u(:,:,jk,jh+nb_ana)
582	ENDDO
583	z2real(:,:)=out_u(:,:,jpk,jh); z2im(:,:)=out_u(:,:,jpk,jh+nb_ana)
584	CALL iom_put( TRIM(tname(jh))//'x_u3d', z3real(:,:,:) )
585	CALL iom_put( TRIM(tname(jh))//'y_u3d', z3im (:,:,:) )
586	CALL iom_put( TRIM(tname(jh))//'x_u2d', z2real(:,:) )
587	CALL iom_put( TRIM(tname(jh))//'y_u2d', z2im (:,:) )
588	z2real(:,:)=out_w(:,:,jpk,jh); z2im(:,:)=out_w(:,:,jpk,jh+nb_ana)
589	CALL iom_put( TRIM(tname(jh))//'x_tabx', z2real(:,:) )
590	CALL iom_put( TRIM(tname(jh))//'y_tabx', z2im (:,:) )
591	# else
592	CALL iom_put( TRIM(tname(jh))//'x_u2d', out_u(:,:,jh) )
593	CALL iom_put( TRIM(tname(jh))//'y_u2d', out_u(:,:,nb_ana+jh) )
594	# endif
595	END DO
596	#endif
597
598	! C) v
599	!/////////
600	#if defined key_dimgout
601	cltext='3d v amplitude and phase; vbar is the last level'
602	CALL dia_wri_dimg(TRIM(cdfile_name_V), TRIM(cltext), out_v, 2*nb_ana, '2')
603	#else
604	# if defined key_3Ddiaharm
605	z3real(:,:,:) = 0._wp; z3im(:,:,:) = 0._wp
606	# endif
607	DO jh = 1, nb_ana
608	# if defined key_3Ddiaharm
609	DO jk=1,jpkm1
610	z3real(:,:,jk)=out_v(:,:,jk,jh)
611	z3im (:,:,jk)=out_v(:,:,jk,jh+nb_ana)
612	ENDDO
613	z2real(:,:)=out_v(:,:,jpk,jh); z2im(:,:)=out_v(:,:,jpk,jh+nb_ana)
614	CALL iom_put( TRIM(tname(jh))//'x_v3d', z3real(:,:,:) )
615	CALL iom_put( TRIM(tname(jh))//'y_v3d', z3im (:,:,:) )
616	CALL iom_put( TRIM(tname(jh))//'x_v2d' , z2real(:,:) )
617	CALL iom_put( TRIM(tname(jh))//'y_v2d' , z2im (:,:) )
618	z2real(:,:)=out_dzi(:,:,jpk,jh); z2im(:,:)=out_dzi(:,:,jpk,jh+nb_ana)
619	CALL iom_put( TRIM(tname(jh))//'x_taby', z2real(:,:) )
620	CALL iom_put( TRIM(tname(jh))//'y_taby', z2im (:,:) )
621	# else
622	CALL iom_put( TRIM(tname(jh))//'x_v2d', out_v(:,:,jh ) )
623	CALL iom_put( TRIM(tname(jh))//'y_v2d', out_v(:,:,jh+nb_ana) )
624	# endif
625	END DO
626
627	#endif
628	! D) w
629	# if defined key_3Ddiaharm
630	# if defined key_dimgout
631	cltext='3d w amplitude and phase; vort_baro is the last level'
632	CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2')
633	# else
634	DO jh = 1, nb_ana
635	DO jk=1,jpkm1
636	z3real(:,:,jk)=out_w(:,:,jk,jh)
637	z3im(:,:,jk)=out_w(:,:,jk,jh+nb_ana)
638	ENDDO
639	CALL iom_put( TRIM(tname(jh))//'x_w3d', z3real(:,:,:) )
640	CALL iom_put( TRIM(tname(jh))//'y_w3d', z3im(:,:,:) )
641	END DO
642	# endif
643
644	! E) dzi + tau_bot
645	# if defined key_dimgout
646	cltext='dzi=g*ro/N2 amplitude and phase'
647	CALL dia_wri_dimg(TRIM(cdfile_name_W), TRIM(cltext), out_w, 2*nb_ana, '2')
648	# else
649	DO jh = 1, nb_ana
650	DO jk=1,jpkm1
651	z3real(:,:,jk)=out_dzi(:,:,jk,jh)
652	z3im(:,:,jk)=out_dzi(:,:,jk,jh+nb_ana)
653	ENDDO
654	CALL iom_put( TRIM(tname(jh))//'x_dzi', z3real(:,:,:) )
655	CALL iom_put( TRIM(tname(jh))//'y_dzi', z3im(:,:,:) )
656	END DO
657	# endif
658	# endif
659
660	!
661	# if defined key_3Ddiaharm
662	DEALLOCATE(z3real, z3im, z2real,z2im)
663	# endif
664
665	END SUBROUTINE dia_wri_harm
666
667
668	SUBROUTINE SUR_DETERMINE_INIT
669	!!---------------------------------------------------------------------------------
670	!! * ROUTINE SUR_DETERMINE_INIT *
671	!!
672	!!---------------------------------------------------------------------------------
673	INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd, jh1_sd, jh2_sd
674	INTEGER :: ipivot
675	REAL(wp) :: zval1, zval2, zcol1, zcol2
676	INTEGER , DIMENSION(jpincomax) :: ipos2
677	!!---------------------------------------------------------------------------------
678	!
679	!
680	ztmp3(:,:) = 0._wp
681	!
682	DO jh1_sd = 1, nsparse
683	DO jh2_sd = 1, nsparse
684	IF( nisparse(jh2_sd) == nisparse(jh1_sd) ) THEN
685	ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) = ztmp3(njsparse(jh1_sd),njsparse(jh2_sd)) &
686	& + valuesparse(jh1_sd)*valuesparse(jh2_sd)
687	ENDIF
688	END DO
689	END DO
690	!
691	DO jj_sd = 1, ninco
692	ipos1(jj_sd) = jj_sd
693	ipos2(jj_sd) = jj_sd
694	END DO
695	!
696	DO ji_sd = 1, ninco
697	!
698	!find greatest non-zero pivot:
699	zval1 = ABS(ztmp3(ji_sd,ji_sd))
700	!
701	ipivot = ji_sd
702	DO jj_sd = ji_sd, ninco
703	zval2 = ABS(ztmp3(ji_sd,jj_sd))
704	IF( zval2 >= zval1 )THEN
705	ipivot = jj_sd
706	zval1 = zval2
707	ENDIF
708	END DO
709	!
710	DO ji1_sd = 1, ninco
711	zcol1 = ztmp3(ji1_sd,ji_sd)
712	zcol2 = ztmp3(ji1_sd,ipivot)
713	ztmp3(ji1_sd,ji_sd) = zcol2
714	ztmp3(ji1_sd,ipivot) = zcol1
715	END DO
716	!
717	ipos2(ji_sd) = ipos1(ipivot)
718	ipos2(ipivot) = ipos1(ji_sd)
719	ipos1(ji_sd) = ipos2(ji_sd)
720	ipos1(ipivot) = ipos2(ipivot)
721	z1_pivot(ji_sd) = 1._wp / ztmp3(ji_sd,ji_sd)
722	DO jj_sd = 1, ninco
723	ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) * z1_pivot(ji_sd)
724	END DO
725	!
726	DO ji2_sd = ji_sd+1, ninco
727	zpilier(ji2_sd,ji_sd) = ztmp3(ji2_sd,ji_sd)
728	DO jj_sd=1,ninco
729	ztmp3(ji2_sd,jj_sd) = ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
730	END DO
731	END DO
732	!
733	END DO
734	!
735	z1_tmp3 = 1._wp / ztmp3(ninco,ninco)
736	!
737	END SUBROUTINE SUR_DETERMINE_INIT
738
739
740	SUBROUTINE SUR_DETERMINE
741	!!---------------------------------------------------------------------------------
742	!! * ROUTINE SUR_DETERMINE *
743	!!
744	!!---------------------------------------------------------------------------------
745	INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd
746	REAL(wp) :: zx1
747	REAL(wp), DIMENSION(jpincomax) :: ztmpx
748	!!---------------------------------------------------------------------------------
749	!
750	DO ji_sd = 1, ninco
751	ztmp4(ji_sd) = ztmp4(ji_sd) * z1_pivot(ji_sd)
752	DO ji2_sd = ji_sd+1, ninco
753	ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
754	END DO
755	END DO
756
757	!system solving:
758	ztmpx(ninco) = ztmp4(ninco) * z1_tmp3
759	DO ji_sd = ninco-1, 1, -1
760	zx1 = 0._wp
761	DO jj_sd = ji_sd+1, ninco
762	zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
763	END DO
764	ztmpx(ji_sd) = ztmp4(ji_sd) - zx1
765	END DO
766
767	DO jj_sd = 1, ninco
768	ztmp7(ipos1(jj_sd)) = ztmpx(jj_sd)
769	END DO
770	!
771	END SUBROUTINE SUR_DETERMINE
772
773
774	!!======================================================================
775	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/r14075_India_uncoupled/src/OCE/DIA/diaharm.F90 @ 15422

Download in other formats: