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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diaharm.F90 in branches/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2015/dev_r5056_CMCC4_simplification/NEMOGCM/NEMO/OPA_SRC/DIA/diaharm.F90 @ 5282

Last change on this file since 5282 was 5282, checked in by diovino, 9 years ago
Dev. branch CMCC4_simplification ticket #1456
File size: 19.2 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 && defined key_tide
9	!!----------------------------------------------------------------------
10	!! 'key_diaharm'
11	!! 'key_tide'
12	!!----------------------------------------------------------------------
13	USE oce ! ocean dynamics and tracers variables
14	USE dom_oce ! ocean space and time domain
15	USE phycst
16	USE dynspg_oce
17	USE dynspg_ts
18	USE daymod
19	USE tide_mod
20	!
21	USE in_out_manager ! I/O units
22	USE iom ! I/0 library
23	USE ioipsl ! NetCDF IPSL library
24	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
25	USE timing ! preformance summary
26	USE wrk_nemo ! working arrays
27
28	IMPLICIT NONE
29	PRIVATE
30
31	LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .TRUE.
32
33	INTEGER, PARAMETER :: jpincomax = 2.*jpmax_harmo
34	INTEGER, PARAMETER :: jpdimsparse = jpincomax30024
35
36	! !! namelist variables
37	INTEGER :: nit000_han ! First time step used for harmonic analysis
38	INTEGER :: nitend_han ! Last time step used for harmonic analysis
39	INTEGER :: nstep_han ! Time step frequency for harmonic analysis
40	INTEGER :: nb_ana ! Number of harmonics to analyse
41
42	INTEGER , ALLOCATABLE, DIMENSION(:) :: name
43	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ana_temp
44	REAL(wp), ALLOCATABLE, DIMENSION(:) :: ana_freq, ut , vt , ft
45	REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: out_eta , out_u, out_v
46
47	INTEGER :: ninco, nsparse
48	INTEGER , DIMENSION(jpdimsparse) :: njsparse, nisparse
49	INTEGER , SAVE, DIMENSION(jpincomax) :: ipos1
50	REAL(wp), DIMENSION(jpdimsparse) :: valuesparse
51	REAL(wp), DIMENSION(jpincomax) :: ztmp4 , ztmp7
52	REAL(wp), SAVE, DIMENSION(jpincomax,jpincomax) :: ztmp3 , zpilier
53	REAL(wp), SAVE, DIMENSION(jpincomax) :: zpivot
54
55	CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: tname ! Names of tidal constituents ('M2', 'K1',...)
56
57	PUBLIC dia_harm ! routine called by step.F90
58
59	!!----------------------------------------------------------------------
60	!! NEMO/OPA 3.5 , NEMO Consortium (2013)
61	!! $Id:$
62	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
63	!!----------------------------------------------------------------------
64	CONTAINS
65
66	SUBROUTINE dia_harm_init
67	!!----------------------------------------------------------------------
68	!! * ROUTINE dia_harm_init *
69	!!
70	!! ** Purpose : Initialization of tidal harmonic analysis
71	!!
72	!! ** Method : Initialize frequency array and nodal factor for nit000_han
73	!!
74	!!--------------------------------------------------------------------
75	INTEGER :: jh, nhan, jk, ji
76	INTEGER :: ios ! Local integer output status for namelist read
77
78	NAMELIST/nam_diaharm/ nit000_han, nitend_han, nstep_han, tname
79	!!----------------------------------------------------------------------
80
81	IF(lwp) THEN
82	WRITE(numout,*)
83	WRITE(numout,*) 'dia_harm_init: Tidal harmonic analysis initialization'
84	WRITE(numout,*) '~~~~~~~ '
85	ENDIF
86	!
87	CALL tide_init_Wave
88	!
89	REWIND( numnam_ref ) ! Namelist nam_diaharm in reference namelist : Tidal harmonic analysis
90	READ ( numnam_ref, nam_diaharm, IOSTAT = ios, ERR = 901)
91	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm in reference namelist', lwp )
92
93	REWIND( numnam_cfg ) ! Namelist nam_diaharm in configuration namelist : Tidal harmonic analysis
94	READ ( numnam_cfg, nam_diaharm, IOSTAT = ios, ERR = 902 )
95	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm in configuration namelist', lwp )
96	IF(lwm) WRITE ( numond, nam_diaharm )
97	!
98	IF(lwp) THEN
99	WRITE(numout,*) 'First time step used for analysis: nit000_han= ', nit000_han
100	WRITE(numout,*) 'Last time step used for analysis: nitend_han= ', nitend_han
101	WRITE(numout,*) 'Time step frequency for harmonic analysis: nstep_han= ', nstep_han
102	ENDIF
103
104	! Basic checks on harmonic analysis time window:
105	! ----------------------------------------------
106	IF( nit000 > nit000_han ) CALL ctl_stop( 'dia_harm_init : nit000_han must be greater than nit000', &
107	& ' restart capability not implemented' )
108	IF( nitend < nitend_han ) CALL ctl_stop( 'dia_harm_init : nitend_han must be lower than nitend', &
109	& 'restart capability not implemented' )
110
111	IF( MOD( nitend_han-nit000_han+1 , nstep_han ) /= 0 ) &
112	& CALL ctl_stop( 'dia_harm_init : analysis time span must be a multiple of nstep_han' )
113
114	nb_ana = 0
115	DO jk=1,jpmax_harmo
116	DO ji=1,jpmax_harmo
117	IF(TRIM(tname(jk)) == Wave(ji)%cname_tide) THEN
118	nb_ana=nb_ana+1
119	ENDIF
120	END DO
121	END DO
122	!
123	IF(lwp) THEN
124	WRITE(numout,*) ' Namelist nam_diaharm'
125	WRITE(numout,*) ' nb_ana = ', nb_ana
126	CALL flush(numout)
127	ENDIF
128	!
129	IF (nb_ana > jpmax_harmo) THEN
130	IF(lwp) WRITE(numout,*) ' E R R O R dia_harm_init : nb_ana must be lower than jpmax_harmo, stop'
131	IF(lwp) WRITE(numout,*) ' jpmax_harmo= ', jpmax_harmo
132	nstop = nstop + 1
133	ENDIF
134
135	ALLOCATE(name (nb_ana))
136	DO jk=1,nb_ana
137	DO ji=1,jpmax_harmo
138	IF (TRIM(tname(jk)) .eq. Wave(ji)%cname_tide) THEN
139	name(jk) = ji
140	EXIT
141	END IF
142	END DO
143	END DO
144
145	! Initialize frequency array:
146	! ---------------------------
147	ALLOCATE( ana_freq(nb_ana), ut(nb_ana), vt(nb_ana), ft(nb_ana) )
148
149	CALL tide_harmo( ana_freq, vt, ut, ft, name, nb_ana )
150
151	IF(lwp) WRITE(numout,*) 'Analysed frequency : ',nb_ana ,'Frequency '
152
153	DO jh = 1, nb_ana
154	IF(lwp) WRITE(numout,*) ' : ',tname(jh),' ',ana_freq(jh)
155	END DO
156
157	! Initialize temporary arrays:
158	! ----------------------------
159	ALLOCATE( ana_temp(jpi,jpj,2*nb_ana,3) )
160	ana_temp(:,:,:,:) = 0._wp
161
162	END SUBROUTINE dia_harm_init
163
164
165	SUBROUTINE dia_harm ( kt )
166	!!----------------------------------------------------------------------
167	!! * ROUTINE dia_harm *
168	!!
169	!! ** Purpose : Tidal harmonic analysis main routine
170	!!
171	!! ** Action : Sums ssh/u/v over time analysis [nit000_han,nitend_han]
172	!!
173	!!--------------------------------------------------------------------
174	INTEGER, INTENT( IN ) :: kt
175	!
176	INTEGER :: ji, jj, jh, jc, nhc
177	REAL(wp) :: ztime, ztemp
178	!!--------------------------------------------------------------------
179	IF( nn_timing == 1 ) CALL timing_start('dia_harm')
180
181	IF( kt == nit000 ) CALL dia_harm_init
182
183	IF( kt >= nit000_han .AND. kt <= nitend_han .AND. MOD(kt,nstep_han) == 0 ) THEN
184
185	ztime = (kt-nit000+1) * rdt
186
187	nhc = 0
188	DO jh = 1, nb_ana
189	DO jc = 1, 2
190	nhc = nhc+1
191	ztemp =( MOD(jc,2) * ft(jh) COS(ana_freq(jh)ztime + vt(jh) + ut(jh)) &
192	& +(1.-MOD(jc,2))* ft(jh) SIN(ana_freq(jh)ztime + vt(jh) + ut(jh)))
193
194	DO jj = 1,jpj
195	DO ji = 1,jpi
196	! Elevation
197	ana_temp(ji,jj,nhc,1) = ana_temp(ji,jj,nhc,1) + ztempsshn(ji,jj) tmask_i(ji,jj)
198	#if defined key_dynspg_ts
199	ana_temp(ji,jj,nhc,2) = ana_temp(ji,jj,nhc,2) + ztempun_b(ji,jj)hur(ji,jj)*umask_i(ji,jj)
200	ana_temp(ji,jj,nhc,3) = ana_temp(ji,jj,nhc,3) + ztempvn_b(ji,jj)hvr(ji,jj)*vmask_i(ji,jj)
201	#endif
202	END DO
203	END DO
204	!
205	END DO
206	END DO
207	!
208	END IF
209
210	IF ( kt == nitend_han ) CALL dia_harm_end
211
212	IF( nn_timing == 1 ) CALL timing_stop('dia_harm')
213
214	END SUBROUTINE dia_harm
215
216
217	SUBROUTINE dia_harm_end
218	!!----------------------------------------------------------------------
219	!! * ROUTINE diaharm_end *
220	!!
221	!! ** Purpose : Compute the Real and Imaginary part of tidal constituents
222	!!
223	!! ** Action : Decompose the signal on the harmonic constituents
224	!!
225	!!--------------------------------------------------------------------
226	INTEGER :: ji, jj, jh, jc, jn, nhan, jl
227	INTEGER :: ksp, kun, keq
228	REAL(wp) :: ztime, ztime_ini, ztime_end
229	REAL(wp) :: X1,X2
230	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ana_amp
231	!!--------------------------------------------------------------------
232	CALL wrk_alloc( jpi , jpj , jpmax_harmo , 2 , ana_amp )
233
234	IF(lwp) WRITE(numout,*)
235	IF(lwp) WRITE(numout,*) 'anharmo_end: kt=nitend_han: Perform harmonic analysis'
236	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
237
238	ztime_ini = nit000_han*rdt ! Initial time in seconds at the beginning of analysis
239	ztime_end = nitend_han*rdt ! Final time in seconds at the end of analysis
240	nhan = (nitend_han-nit000_han+1)/nstep_han ! Number of dumps used for analysis
241
242	ninco = 2*nb_ana
243
244	ksp = 0
245	keq = 0
246	DO jn = 1, nhan
247	ztime=( (nhan-jn)ztime_ini + (jn-1)ztime_end )/FLOAT(nhan-1)
248	keq = keq + 1
249	kun = 0
250	DO jh = 1, nb_ana
251	DO jc = 1, 2
252	kun = kun + 1
253	ksp = ksp + 1
254	nisparse(ksp) = keq
255	njsparse(ksp) = kun
256	valuesparse(ksp) = ( MOD(jc,2) * ft(jh) * COS(ana_freq(jh)*ztime + vt(jh) + ut(jh)) &
257	& + (1.-MOD(jc,2))* ft(jh) * SIN(ana_freq(jh)*ztime + vt(jh) + ut(jh)) )
258	END DO
259	END DO
260	END DO
261
262	nsparse = ksp
263
264	! Elevation:
265	DO jj = 1, jpj
266	DO ji = 1, jpi
267	! Fill input array
268	kun = 0
269	DO jh = 1, nb_ana
270	DO jc = 1, 2
271	kun = kun + 1
272	ztmp4(kun)=ana_temp(ji,jj,kun,1)
273	END DO
274	END DO
275
276	CALL SUR_DETERMINE(jj)
277
278	! Fill output array
279	DO jh = 1, nb_ana
280	ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
281	ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
282	END DO
283	END DO
284	END DO
285
286	ALLOCATE( out_eta(jpi,jpj,2*nb_ana), &
287	& out_u (jpi,jpj,2*nb_ana), &
288	& out_v (jpi,jpj,2*nb_ana) )
289
290	DO jj = 1, jpj
291	DO ji = 1, jpi
292	DO jh = 1, nb_ana
293	X1 = ana_amp(ji,jj,jh,1)
294	X2 =-ana_amp(ji,jj,jh,2)
295	out_eta(ji,jj,jh ) = X1 * tmask_i(ji,jj)
296	out_eta(ji,jj,jh+nb_ana) = X2 * tmask_i(ji,jj)
297	END DO
298	END DO
299	END DO
300
301	! ubar:
302	DO jj = 1, jpj
303	DO ji = 1, jpi
304	! Fill input array
305	kun=0
306	DO jh = 1,nb_ana
307	DO jc = 1,2
308	kun = kun + 1
309	ztmp4(kun)=ana_temp(ji,jj,kun,2)
310	END DO
311	END DO
312
313	CALL SUR_DETERMINE(jj+1)
314
315	! Fill output array
316	DO jh = 1, nb_ana
317	ana_amp(ji,jj,jh,1) = ztmp7((jh-1)*2+1)
318	ana_amp(ji,jj,jh,2) = ztmp7((jh-1)*2+2)
319	END DO
320
321	END DO
322	END DO
323
324	DO jj = 1, jpj
325	DO ji = 1, jpi
326	DO jh = 1, nb_ana
327	X1= ana_amp(ji,jj,jh,1)
328	X2=-ana_amp(ji,jj,jh,2)
329	out_u(ji,jj, jh) = X1 * umask_i(ji,jj)
330	out_u(ji,jj,nb_ana+jh) = X2 * umask_i(ji,jj)
331	ENDDO
332	ENDDO
333	ENDDO
334
335	! vbar:
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	ztmp4(kun)=ana_temp(ji,jj,kun,3)
344	END DO
345	END DO
346
347	CALL SUR_DETERMINE(jj+1)
348
349	! Fill output array
350	DO jh = 1, nb_ana
351	ana_amp(ji,jj,jh,1)=ztmp7((jh-1)*2+1)
352	ana_amp(ji,jj,jh,2)=ztmp7((jh-1)*2+2)
353	END DO
354
355	END DO
356	END DO
357
358	DO jj = 1, jpj
359	DO ji = 1, jpi
360	DO jh = 1, nb_ana
361	X1=ana_amp(ji,jj,jh,1)
362	X2=-ana_amp(ji,jj,jh,2)
363	out_v(ji,jj, jh)=X1 * vmask_i(ji,jj)
364	out_v(ji,jj,nb_ana+jh)=X2 * vmask_i(ji,jj)
365	END DO
366	END DO
367	END DO
368
369	CALL dia_wri_harm ! Write results in files
370	CALL wrk_dealloc( jpi , jpj , jpmax_harmo , 2 , ana_amp )
371	!
372	END SUBROUTINE dia_harm_end
373
374
375	SUBROUTINE dia_wri_harm
376	!!--------------------------------------------------------------------
377	!! * ROUTINE dia_wri_harm *
378	!!
379	!! ** Purpose : Write tidal harmonic analysis results in a netcdf file
380	!!--------------------------------------------------------------------
381	CHARACTER(LEN=lc) :: cltext
382	CHARACTER(LEN=lc) :: &
383	cdfile_name_T , & ! name of the file created (T-points)
384	cdfile_name_U , & ! name of the file created (U-points)
385	cdfile_name_V ! name of the file created (V-points)
386	INTEGER :: jh
387	!!----------------------------------------------------------------------
388
389	!#if
390	!#endif
391
392	IF(lwp) WRITE(numout,*) ' '
393	IF(lwp) WRITE(numout,*) 'dia_wri_harm : Write harmonic analysis results'
394	!#if
395	!#endif
396	IF(lwp) WRITE(numout,*) ' '
397
398	! A) Elevation
399	!/////////////
400	!
401	!#if
402	!#else
403	DO jh = 1, nb_ana
404	CALL iom_put( TRIM(tname(jh))//'x', out_eta(:,:,jh) )
405	CALL iom_put( TRIM(tname(jh))//'y', out_eta(:,:,nb_ana+jh) )
406	END DO
407	!#endif
408
409	! B) ubar
410	!/////////
411	!
412	!#if
413	!#else
414	DO jh = 1, nb_ana
415	CALL iom_put( TRIM(tname(jh))//'x_u', out_u(:,:,jh) )
416	CALL iom_put( TRIM(tname(jh))//'y_u', out_u(:,:,nb_ana+jh) )
417	END DO
418	!#endif
419
420	! C) vbar
421	!/////////
422	!
423	!#if
424	!#else
425	DO jh = 1, nb_ana
426	CALL iom_put( TRIM(tname(jh))//'x_v', out_v(:,:,jh ) )
427	CALL iom_put( TRIM(tname(jh))//'y_v', out_v(:,:,jh+nb_ana) )
428	END DO
429	!#endif
430	!
431	END SUBROUTINE dia_wri_harm
432
433
434	SUBROUTINE SUR_DETERMINE(init)
435	!!---------------------------------------------------------------------------------
436	!! * ROUTINE SUR_DETERMINE *
437	!!
438	!!
439	!!
440	!!---------------------------------------------------------------------------------
441	INTEGER, INTENT(in) :: init
442	!
443	INTEGER :: ji_sd, jj_sd, ji1_sd, ji2_sd, jk1_sd, jk2_sd
444	REAL(wp) :: zval1, zval2, zx1
445	REAL(wp), POINTER, DIMENSION(:) :: ztmpx, zcol1, zcol2
446	INTEGER , POINTER, DIMENSION(:) :: ipos2, ipivot
447	!---------------------------------------------------------------------------------
448	CALL wrk_alloc( jpincomax , ztmpx , zcol1 , zcol2 )
449	CALL wrk_alloc( jpincomax , ipos2 , ipivot )
450
451	IF( init == 1 ) THEN
452	IF( nsparse > jpdimsparse ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : nsparse .GT. jpdimsparse')
453	IF( ninco > jpincomax ) CALL ctl_stop( 'STOP', 'SUR_DETERMINE : ninco .GT. jpincomax')
454	!
455	ztmp3(:,:) = 0._wp
456	!
457	DO jk1_sd = 1, nsparse
458	DO jk2_sd = 1, nsparse
459	nisparse(jk2_sd) = nisparse(jk2_sd)
460	njsparse(jk2_sd) = njsparse(jk2_sd)
461	IF( nisparse(jk2_sd) == nisparse(jk1_sd) ) THEN
462	ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) = ztmp3(njsparse(jk1_sd),njsparse(jk2_sd)) &
463	& + valuesparse(jk1_sd)*valuesparse(jk2_sd)
464	ENDIF
465	END DO
466	END DO
467	!
468	DO jj_sd = 1 ,ninco
469	ipos1(jj_sd) = jj_sd
470	ipos2(jj_sd) = jj_sd
471	END DO
472	!
473	DO ji_sd = 1 , ninco
474	!
475	!find greatest non-zero pivot:
476	zval1 = ABS(ztmp3(ji_sd,ji_sd))
477	!
478	ipivot(ji_sd) = ji_sd
479	DO jj_sd = ji_sd, ninco
480	zval2 = ABS(ztmp3(ji_sd,jj_sd))
481	IF( zval2.GE.zval1 )THEN
482	ipivot(ji_sd) = jj_sd
483	zval1 = zval2
484	ENDIF
485	END DO
486	!
487	DO ji1_sd = 1, ninco
488	zcol1(ji1_sd) = ztmp3(ji1_sd,ji_sd)
489	zcol2(ji1_sd) = ztmp3(ji1_sd,ipivot(ji_sd))
490	ztmp3(ji1_sd,ji_sd) = zcol2(ji1_sd)
491	ztmp3(ji1_sd,ipivot(ji_sd)) = zcol1(ji1_sd)
492	END DO
493	!
494	ipos2(ji_sd) = ipos1(ipivot(ji_sd))
495	ipos2(ipivot(ji_sd)) = ipos1(ji_sd)
496	ipos1(ji_sd) = ipos2(ji_sd)
497	ipos1(ipivot(ji_sd)) = ipos2(ipivot(ji_sd))
498	zpivot(ji_sd) = ztmp3(ji_sd,ji_sd)
499	DO jj_sd = 1, ninco
500	ztmp3(ji_sd,jj_sd) = ztmp3(ji_sd,jj_sd) / zpivot(ji_sd)
501	END DO
502	!
503	DO ji2_sd = ji_sd+1, ninco
504	zpilier(ji2_sd,ji_sd)=ztmp3(ji2_sd,ji_sd)
505	DO jj_sd=1,ninco
506	ztmp3(ji2_sd,jj_sd)= ztmp3(ji2_sd,jj_sd) - ztmp3(ji_sd,jj_sd) * zpilier(ji2_sd,ji_sd)
507	END DO
508	END DO
509	!
510	END DO
511	!
512	ENDIF ! End init==1
513
514	DO ji_sd = 1, ninco
515	ztmp4(ji_sd) = ztmp4(ji_sd) / zpivot(ji_sd)
516	DO ji2_sd = ji_sd+1, ninco
517	ztmp4(ji2_sd) = ztmp4(ji2_sd) - ztmp4(ji_sd) * zpilier(ji2_sd,ji_sd)
518	END DO
519	END DO
520
521	!system solving:
522	ztmpx(ninco) = ztmp4(ninco) / ztmp3(ninco,ninco)
523	ji_sd = ninco
524	DO ji_sd = ninco-1, 1, -1
525	zx1 = 0._wp
526	DO jj_sd = ji_sd+1, ninco
527	zx1 = zx1 + ztmpx(jj_sd) * ztmp3(ji_sd,jj_sd)
528	END DO
529	ztmpx(ji_sd) = ztmp4(ji_sd)-zx1
530	END DO
531
532	DO jj_sd =1, ninco
533	ztmp7(ipos1(jj_sd))=ztmpx(jj_sd)
534	END DO
535
536	CALL wrk_dealloc( jpincomax , ztmpx , zcol1 , zcol2 )
537	CALL wrk_dealloc( jpincomax , ipos2 , ipivot )
538	!
539	END SUBROUTINE SUR_DETERMINE
540
541	#else
542	!!----------------------------------------------------------------------
543	!! Default case : Empty module
544	!!----------------------------------------------------------------------
545	LOGICAL, PUBLIC, PARAMETER :: lk_diaharm = .FALSE.
546	CONTAINS
547	SUBROUTINE dia_harm ( kt ) ! Empty routine
548	INTEGER, INTENT( IN ) :: kt
549	WRITE(,) 'dia_harm: you should not have seen this print'
550	END SUBROUTINE dia_harm
551	#endif
552
553	!!======================================================================
554	END MODULE diaharm

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

Download in other formats: