Context Navigation

diaharm_fast.F90

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

Line
1	MODULE diaharm_fast
2	!!======================================================================
3	!! * MODULE example *
4	!! Ocean physics: On line harmonic analyser
5	!!
6	!!=====================================================================
7
8	#if defined key_diaharm_fast
9
10	!!----------------------------------------------------------------------
11	!! 'key_harm_ana' : Calculate harmonic analysis
12	!!----------------------------------------------------------------------
13	!! harm_ana :
14	!! harm_ana_init :
15	!! NB: 2017-12 : add 3D harmonic analysis of velocities
16	!! integration of Maria Luneva's development
17	!! 'key_3Ddiaharm'
18	!!----------------------------------------------------------------------
19
20	USE oce ! ocean dynamics and tracers
21	USE dom_oce ! ocean space and time domain
22	USE iom
23	USE in_out_manager ! I/O units
24	USE phycst ! physical constants
25	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
26	USE bdy_oce ! ocean open boundary conditions
27	USE bdytides ! tidal bdy forcing
28	USE daymod ! calendar
29	USE tideini
30	USE restart
31	USE ioipsl, ONLY : ju2ymds ! for calendar
32	!
33	!
34	USE timing ! preformance summary
35	USE zdf_oce
36
37	IMPLICIT NONE
38	PRIVATE
39
40	!! * Routine accessibility
41	PUBLIC dia_harm_fast ! routine called in step.F90 module
42	LOGICAL, PUBLIC, PARAMETER :: lk_diaharm_fast = .TRUE. ! to be run or not
43	LOGICAL, PUBLIC :: lk_diaharm_2D ! = .TRUE. ! to run 2d
44	LOGICAL, PUBLIC :: lk_diaharm_3D ! = .TRUE. ! to run 3d
45
46	!! * Module variables
47	INTEGER, PARAMETER :: nharm_max = jpmax_harmo ! max number of harmonics to be analysed
48	INTEGER, PARAMETER :: nhm_max = 2*nharm_max+1
49	INTEGER, PARAMETER :: nvab = 2 ! number of 3D variables
50	INTEGER :: nharm
51	INTEGER :: nhm
52	INTEGER :: & !!! toto namelist (namtoto)
53	nflag = 1 ! default value of nflag
54	REAL(wp), DIMENSION(nharm_max) :: &
55	om_tide ! tidal frequencies ( rads/sec)
56	REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:) :: &
57	bzz,c,x ! work arrays
58	REAL(wp) :: cca,ssa,zm,bt,dd_cumul
59	!
60	REAL(wp), PUBLIC :: fjulday_startharm !: Julian Day since start of harmonic analysis
61	REAL(wp), PUBLIC, ALLOCATABLE,DIMENSION(:) :: anau, anav, anaf ! nodel/phase corrections used by diaharmana
62	REAL(WP), ALLOCATABLE,SAVE,DIMENSION(:,:) :: cc,a
63	!
64	INTEGER :: nvar_2d, nvar_3d !: number of 2d and 3d variables to analyse
65	INTEGER, ALLOCATABLE,DIMENSION(:) :: m_posi_2d, m_posi_3d
66
67	! Name of variables used in the restart
68	CHARACTER( LEN = 10 ), DIMENSION(5), PARAMETER :: m_varName2d = (/'ssh','u2d','v2d','ubfr','vbfr'/)
69	CHARACTER( LEN = 10 ), DIMENSION(4), PARAMETER :: m_varName3d = (/'rho','u3d','v3d','w3d'/)
70	!
71	REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,: ) :: g_cosamp2D, g_sinamp2D, g_cumul_var2D
72	REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:,:,:) :: g_cosamp3D, g_sinamp3D, g_cumul_var3D
73	!
74	REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:) :: g_out2D,h_out2D ! arrays for output
75	REAL(wp), ALLOCATABLE,SAVE,DIMENSION(:,:,:) :: g_out3D,h_out3D ! arrays for 3D output
76	!
77	! NAMELIST
78	LOGICAL, PUBLIC :: ln_diaharm_store !: =T Stores data for harmonic Analysis
79	LOGICAL, PUBLIC :: ln_diaharm_compute !: =T Compute harmonic Analysis
80	LOGICAL, PUBLIC :: ln_diaharm_read_restart !: =T Read restart from a previous run
81	LOGICAL, PUBLIC :: ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d
82	INTEGER :: nb_ana ! Number of harmonics to analyse
83	CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: tname ! Names of tidal constituents ('M2', 'K1',...)
84	INTEGER , ALLOCATABLE, DIMENSION(:) :: ntide_all ! INDEX within the full set of constituents (tide.h90)
85	INTEGER , ALLOCATABLE, DIMENSION(:) :: ntide_sub ! INDEX within the subset of constituents pass in input
86
87	!! * Substitutions
88
89	!!----------------------------------------------------------------------
90	!! OPA 9.0 , LOCEAN-IPSL (2005)
91	!! or LIM 2.0 , UCL-LOCEAN-IPSL (2005)
92	!! or TOP 1.0 , LOCEAN-IPSL (2005)
93	!! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/module_example,v 1.3 2005/03/27 18:34:47 opalod Exp $
94	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
95	!!----------------------------------------------------------------------
96
97	CONTAINS
98
99	SUBROUTINE dia_harm_fast( kt )
100	!!----------------------------------------------------------------------
101	!! * ROUTINE harm_ana *
102	!!
103	!! ** Purpose : Harmonic analyser
104	!!
105	!! ** Method :
106	!!
107	!! ** Action : - first action (share memory array/varible modified
108	!! in this routine
109	!! - second action .....
110	!! - .....
111	!!
112	!! References :
113	!! Give references if exist otherwise suppress these lines
114	!!
115	!! History :
116	!! 9.0 ! 03-08 (Autor Names) Original code
117	!! ! 02-08 (Author names) brief description of modifications
118	!!----------------------------------------------------------------------
119	!! * Modules used
120
121	!! * arguments
122	INTEGER, INTENT( in ) :: &
123	kt ! describe it!!!
124
125	!! * local declarations
126	INTEGER :: ji, jk, jj ! dummy loop arguments
127	INTEGER :: jh, i1, i2, jgrid
128	INTEGER :: j2d, j3d
129	REAL(WP) :: sec2start
130	!!--------------------------------------------------------------------
131
132	IF( nn_timing == 1 ) CALL timing_start( 'dia_harm_fast' )
133	IF( kt == nit000 ) CALL harm_ana_init ! Initialization (first time-step only)
134
135	IF ( ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
136
137	! this bit done every time step
138	nhm=2*nb_ana+1
139	c(1) = 1.0
140
141	sec2start = nint( (fjulday-fjulday_startharm)*86400._wp )
142
143	DO jh=1,nb_ana
144	c(2jh ) = anaf(jh)cos( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
145	c(2jh+1) = anaf(jh)sin( sec2start*om_tide(jh) + anau(jh) + anav(jh) )
146	ENDDO
147
148	! CUMULATE
149	DO ji=1,jpi ! loop lon
150	DO jj=1,jpj ! loop lat
151	DO jh=1,nhm ! loop harmonic
152
153	DO j2d=1,nvar_2d
154	IF ( m_posi_2d(j2d) .eq. 1 ) dd_cumul = c(jh) * sshn(ji,jj) * ssmask (ji,jj) ! analysis elevation
155	IF ( m_posi_2d(j2d) .eq. 2 ) dd_cumul = c(jh) * un_b(ji,jj) * ssumask(ji,jj) ! analysis depth average velocities
156	IF ( m_posi_2d(j2d) .eq. 3 ) dd_cumul = c(jh) * vn_b(ji,jj) * ssvmask(ji,jj)
157	IF ( m_posi_2d(j2d) .eq. 4 ) dd_cumul = c(jh) * bfrua(ji,jj) * un(ji,jj,mbku(ji,jj)) * ssumask(ji,jj) ! analysis bottom friction
158	IF ( m_posi_2d(j2d) .eq. 5 ) dd_cumul = c(jh) * bfrva(ji,jj) * vn(ji,jj,mbkv(ji,jj)) * ssvmask(ji,jj)
159	g_cumul_var2D(jh,ji,jj,j2d) = g_cumul_var2D(jh,ji,jj,j2d) + dd_cumul
160	ENDDO
161
162	DO j3d=1,nvar_3d
163	DO jk=1,jpkm1
164	IF ( m_posi_3d(j3d) .eq. 1 ) dd_cumul = c(jh) * rhd(ji,jj,jk) * tmask(ji,jj,jk)
165	IF ( m_posi_3d(j3d) .eq. 2 ) dd_cumul = c(jh) * ( un(ji,jj,jk)-un_b(ji,jj) ) * umask(ji,jj,jk)
166	IF ( m_posi_3d(j3d) .eq. 3 ) dd_cumul = c(jh) * ( vn(ji,jj,jk)-vn_b(ji,jj) ) * vmask(ji,jj,jk)
167	IF ( m_posi_3d(j3d) .eq. 4 ) dd_cumul = c(jh) * wn(ji,jj,jk) * wmask(ji,jj,jk)
168	g_cumul_var3D(jh,ji,jj,jk,j3d) = g_cumul_var3D(jh,ji,jj,jk,j3d) + dd_cumul
169	ENDDO
170	ENDDO
171
172	ENDDO ! end loop harmonic
173	ENDDO ! end loop lat
174	ENDDO ! end loop lon
175
176	! Compute nodal factor cumulative cross-product
177	DO i1=1,nhm
178	DO i2=1,nhm
179	cc(i1,i2)=cc(i1,i2)+c(i1)*c(i2)
180	ENDDO
181	ENDDO
182
183	! Output RESTART
184	IF( kt == nitrst ) THEN
185	CALL harm_rst_write(kt) ! Dump out data for a restarted run
186	ENDIF
187
188	! At End of run
189	IF ( kt == nitend ) THEN
190
191	IF(lwp) WRITE(numout,*)
192	IF(lwp) WRITE(numout,*) 'harm_ana : harmonic analysis of tides at end of run'
193	IF(lwp) WRITE(numout,*) '~~~~~~~~~'
194
195	IF( ln_diaharm_compute ) THEN
196
197	! INITIALISE TABLE TO 0
198	IF ( nvar_2d .gt. 0 ) THEN
199	g_cosamp2D = 0.0_wp
200	g_sinamp2D = 0.0_wp
201	ENDIF
202	IF ( nvar_3d .gt. 0 ) THEN
203	g_cosamp3D = 0.0_wp
204	g_sinamp3D = 0.0_wp
205	ENDIF
206
207	! FIRST OUTPUT 2D VARIABLES
208	DO jgrid=1,nvar_2d ! loop number of 2d variables (ssh, U2d, V2d, UVfric) to analyse harmonically
209	DO ji=1,jpi ! loop lon
210	DO jj=1,jpj ! loop lat
211	bt = 1.0_wp; bzz(:) = 0.0_wp
212	DO jh=1,nhm ! loop harmonic
213	bzz(jh) = g_cumul_var2D(jh,ji,jj,jgrid)
214	bt = bt*bzz(jh)
215	ENDDO
216	! Copy back original cumulated nodal factor
217	a(:,:) = cc(:,:)
218	! now do gaussian elimination of the system
219	! a * x = b
220	! the matrix x is (a0,a1,b1,a2,b2 ...)
221	! the matrix a and rhs b solved here for x
222	x=0.0_wp
223	IF(bt.ne.0.) THEN
224	CALL gelim( a, bzz, x, nhm )
225	! Backup output in variables
226	DO jh=1,nb_ana
227	g_cosamp2D(jh,ji,jj,jgrid) = x(jh*2 )
228	g_sinamp2D(jh,ji,jj,jgrid) = x(jh*2+1)
229	ENDDO
230	g_cosamp2D( 0,ji,jj,jgrid) = x(1)
231	g_sinamp2D( 0,ji,jj,jgrid) = 0.0_wp
232	ENDIF ! bt.ne.0.
233	ENDDO ! jj
234	ENDDO ! ji
235	ENDDO ! jgrid
236
237	! SECOND OUTPUT 3D VARIABLES
238	DO jgrid=1,nvar_3d ! loop number of 3d variables rho, U, V, W
239	DO jk=1,jpkm1 ! loop over vertical level
240	DO ji=1,jpi ! loop over lon
241	DO jj=1,jpj ! loop over lat
242	bt = 1.0_wp; bzz(:) = 0.0_wp
243	DO jh=1,nhm
244	bzz(jh) = g_cumul_var3D(jh,ji,jj,jk,jgrid)
245	bt = bt*bzz(jh)
246	ENDDO
247	! Copy back original cumulated nodal factor
248	a(:,:) = cc(:,:)
249	! now do gaussian elimination of the system
250	! a * x = b
251	! the matrix x is (a0,a1,b1,a2,b2 ...)
252	! the matrix a and rhs b solved here for x
253	x=0.0_wp
254	IF(bt.ne.0.) THEN
255	CALL gelim( a, bzz, x, nhm )
256	! Backup output in variables
257	DO jh=1,nb_ana
258	g_cosamp3D(jh,ji,jj,jk,jgrid) = x(jh*2 )
259	g_sinamp3D(jh,ji,jj,jk,jgrid) = x(jh*2+1)
260	ENDDO
261	g_cosamp3D ( 0,ji,jj,jk,jgrid) = x(1)
262	g_sinamp3D ( 0,ji,jj,jk,jgrid) = 0.0_wp
263	ENDIF ! bt.ne.0.
264	ENDDO ! jj
265	ENDDO ! ji
266	ENDDO ! jk
267	ENDDO ! jgrid
268
269	CALL harm_ana_out ! output analysis (last time step)
270
271	ELSE ! ln_harmana_compute = False
272	IF(lwp) WRITE(numout,*) " Skipping Computing harmonics at last step"
273
274	ENDIF ! ln_harmana_compute
275	ENDIF ! kt == nitend
276
277	ENDIF
278
279	IF( nn_timing == 1 ) CALL timing_stop( 'dia_harm_fast' )
280
281	END SUBROUTINE dia_harm_fast
282
283	SUBROUTINE harm_ana_init
284	!!----------------------------------------------------------------------
285	!! * ROUTINE harm_ana_init *
286	!!
287	!! ** Purpose : initialization of ....
288	!!
289	!! ** Method : blah blah blah ...
290	!!
291	!! ** input : Namlist namexa
292	!!
293	!! ** Action : ...
294	!!
295	!! history :
296	!! 9.0 ! 03-08 (Autor Names) Original code
297	!!----------------------------------------------------------------------
298	!! * local declarations
299	INTEGER :: ji, jk, jh ! dummy loop indices
300	INTEGER :: ios ! Local integer output status for namelist read
301	INTEGER :: k2d, k3d ! dummy number of analysis
302	NAMELIST/nam_diaharm_fast/ ln_diaharm_store, ln_diaharm_compute, ln_diaharm_read_restart, ln_ana_ssh, ln_ana_uvbar, ln_ana_bfric, ln_ana_rho, ln_ana_uv3d, ln_ana_w3d, tname
303	!!----------------------------------------------------------------------
304
305	lk_diaharm_2D = .TRUE. ! to run 2d
306	lk_diaharm_3D = .TRUE. ! to run 3d
307
308	IF(lwp) WRITE(numout,*)
309	IF(lwp) WRITE(numout,*) 'harm_init : initialization of harmonic analysis of tides'
310	IF(lwp) WRITE(numout,*) '~~~~~~~~~'
311
312	! GET NAMELIST DETAILS
313	REWIND( numnam_ref ) ! Namelist nam_diaharm_fast in reference namelist : Tidal harmonic analysis
314	READ ( numnam_ref, nam_diaharm_fast, IOSTAT = ios, ERR = 901)
315	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in reference namelist', lwp )
316
317	REWIND( numnam_cfg ) ! Namelist nam_diaharm_fast in configuration namelist : Tidal harmonic analysis
318	READ ( numnam_cfg, nam_diaharm_fast, IOSTAT = ios, ERR = 902 )
319	902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_diaharm_fast in configuration namelist', lwp )
320	IF(lwm) WRITE ( numond, nam_diaharm_fast )
321
322	! GET NUMBER OF HARMONIC TO ANALYSE - from diaharm.F90
323	nb_ana = 0
324	DO jk=1,jpmax_harmo
325	DO ji=1,nb_harmo
326	IF(TRIM(tname(jk)) == Wave( ntide(ji) )%cname_tide ) THEN
327	nb_ana=nb_ana+1
328	ENDIF
329	END DO
330	END DO
331	!
332	IF(lwp) THEN
333	WRITE(numout,*) ' Namelist nam_diaharm_fast'
334	WRITE(numout,*) ' nb_ana = ', nb_ana
335	CALL flush(numout)
336	ENDIF
337	!
338	IF (nb_ana > nharm_max) THEN
339	IF(lwp) WRITE(numout,*) ' E R R O R harm_ana : nb_ana must be lower than nharm_max, stop'
340	IF(lwp) WRITE(numout,*) ' nharm_max = ', nharm_max
341	nstop = nstop + 1
342	ENDIF
343
344	ALLOCATE(ntide_all(nb_ana))
345	ALLOCATE(ntide_sub(nb_ana))
346
347	DO jk=1,nb_ana
348	DO ji=1,nb_harmo
349	IF (TRIM(tname(jk)) .eq. Wave( ntide(ji) )%cname_tide ) THEN
350	ntide_sub(jk) = ji
351	ntide_all(jk) = ntide(ji)
352	EXIT
353	END IF
354	END DO
355	END DO
356
357	! SEARCH HOW MANY VARIABLES 2D AND 3D TO PROCESS
358	nvar_2d = 0; nvar_3d = 0
359	IF ( ln_ana_ssh ) nvar_2d = nvar_2d + 1 ! analysis elevation
360	IF ( ln_ana_uvbar ) nvar_2d = nvar_2d + 2 ! analysis depth-averaged velocity
361	IF ( ln_ana_bfric ) nvar_2d = nvar_2d + 2 ! analysis bottom friction
362
363	IF ( ln_ana_rho ) nvar_3d = nvar_3d + 1 ! analysis density
364	IF ( ln_ana_uv3d ) nvar_3d = nvar_3d + 2 ! analysis 3D horizontal velocities
365	IF ( ln_ana_w3d ) nvar_3d = nvar_3d + 1 ! analysis 3D vertical velocity
366
367	! CHECK IF SOMETHING TO RUN
368	IF ( nvar_2d .eq. 0 ) lk_diaharm_2D = .FALSE. ! no 2d to run
369	IF ( nvar_3d .eq. 0 ) lk_diaharm_3D = .FALSE. ! no 3d to run
370
371	IF ( ln_diaharm_store .and. ( lk_diaharm_2D .or. lk_diaharm_3D) ) THEN
372
373	! DO ALLOCATIONS
374	IF ( lk_diaharm_2D ) THEN
375	ALLOCATE( g_cumul_var2D(nb_ana*2+1,jpi,jpj, nvar_2d) )
376	ALLOCATE( g_cosamp2D( 0:nb_ana*2+1,jpi,jpj, nvar_2d) )
377	ALLOCATE( g_sinamp2D( 0:nb_ana*2+1,jpi,jpj, nvar_2d) )
378	ALLOCATE( g_out2D (jpi,jpj) )
379	ALLOCATE( h_out2D (jpi,jpj) )
380	ALLOCATE( m_posi_2d( nvar_2d ) ); m_posi_2d(:)=0
381	ENDIF
382
383	IF ( lk_diaharm_3D ) THEN
384	ALLOCATE( g_cumul_var3D(nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
385	ALLOCATE( g_cosamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
386	ALLOCATE( g_sinamp3D( 0:nb_ana*2+1,jpi,jpj,jpk,nvar_3d) )
387	ALLOCATE( g_out3D (jpi,jpj,jpk) )
388	ALLOCATE( h_out3D (jpi,jpj,jpk) )
389	ALLOCATE( m_posi_3d( nvar_3d ) ); m_posi_3d(:)=0
390	ENDIF
391
392	ALLOCATE( cc(nb_ana2+1,nb_ana2+1) )
393	ALLOCATE( a (nb_ana2+1,nb_ana2+1) )
394	ALLOCATE( bzz(nb_ana*2+1) )
395	ALLOCATE( x (nb_ana*2+1) )
396	ALLOCATE( c (nb_ana*2+1) )
397	ALLOCATE( anau(nb_ana) )
398	ALLOCATE( anav(nb_ana) )
399	ALLOCATE( anaf(nb_ana) )
400	! END ALLOCATE
401
402	! STORE INDEX OF WHAT TO PRODUCE DEPENDING ON ACTIVATED LOGICAL
403	! MAKES THINGS EASIER AND FASTER LATER
404	! !!! UGLY !!!
405	jh = 1; k2d = 0;
406	IF ( ln_ana_ssh ) THEN
407	k2d = k2d + 1; m_posi_2d(k2d) = jh
408	IF(lwp) WRITE(numout,*) " - ssh harmonic analysis activated (ln_ana_ssh)"
409	ENDIF
410	jh = jh + 1
411	IF ( ln_ana_uvbar ) THEN
412	k2d = k2d + 1; m_posi_2d(k2d) = jh
413	jh = jh + 1
414	k2d = k2d + 1; m_posi_2d(k2d) = jh
415	IF(lwp) WRITE(numout,*) " - barotropic currents harmonic analysis activated (ln_ana_uvbar)"
416	ELSE
417	jh = jh + 1
418	ENDIF
419	jh = jh + 1
420	IF ( ln_ana_bfric ) THEN
421	k2d = k2d + 1; m_posi_2d(k2d) = jh
422	jh = jh + 1;
423	k2d = k2d + 1; m_posi_2d(k2d) = jh
424	IF(lwp) WRITE(numout,*) " - bottom friction harmonic analysis activated (ln_ana_vbfr)"
425	ELSE
426	jh = jh + 1
427	ENDIF
428
429	! and for 3D
430	jh = 1; k3d = 0;
431	IF ( ln_ana_rho ) THEN
432	k3d = k3d + 1; m_posi_3d(k3d) = jh
433	IF(lwp) WRITE(numout,*) " - 3D density harmonic analysis activated (ln_ana_rho)"
434	ENDIF
435	jh = jh + 1
436	IF ( ln_ana_uv3d ) THEN
437	k3d = k3d + 1; m_posi_3d(k3d) = jh
438	jh = jh + 1
439	k3d = k3d + 1; m_posi_3d(k3d) = jh
440	IF(lwp) WRITE(numout,*) " - 3D horizontal currents harmonic analysis activated (ln_ana_uv3d)"
441	ELSE
442	jh = jh + 1
443	ENDIF
444	jh = jh + 1
445	IF ( ln_ana_w3d ) THEN
446	k3d = k3d + 1; m_posi_3d(k3d) = jh
447	IF(lwp) WRITE(numout,*) " - 3D vertical currents harmonic analysis activated (ln_ana_w3d)"
448	ENDIF
449
450	! SELECT AND STORE FREQUENCIES
451	IF(lwp) WRITE(numout,*) 'Analysed frequency : ',nb_ana ,'Frequency '
452	DO jh=1,nb_ana
453	om_tide(jh) = omega_tide( ntide_sub(jh) )
454	IF(lwp) WRITE(numout,*) ' - ',tname(jh),' ',om_tide(jh)
455	ENDDO
456
457	! READ RESTART IF
458	IF ( ln_diaharm_read_restart ) THEN
459	IF (lwp) WRITE(numout,*) "Reading previous harmonic data from previous run"
460	! Need to read in bssh bz, cc anau anav and anaf
461	call harm_rst_read ! This reads in from the previous day
462	! Currrently the data in in assci format
463	ELSE
464
465	IF (lwp) WRITE(numout,*) "Starting harmonic analysis from Fresh "
466
467	IF ( lk_diaharm_2D ) g_cumul_var2D(:,:,:,: ) = 0.0_wp
468	IF ( lk_diaharm_3D ) g_cumul_var3D(:,:,:,:,:) = 0.0_wp
469	cc = 0.0_wp
470	a (:,:) = 0.0_wp ! NB
471	bzz (:) = 0.0_wp
472	x (:) = 0.0_wp
473	c (:) = 0.0_wp
474	anau (:) = 0.0_wp
475	anav (:) = 0.0_wp
476	anaf (:) = 0.0_wp
477
478	DO jh = 1, nb_ana
479	anau(jh) = utide ( ntide_sub(jh) )
480	anav(jh) = v0tide( ntide_sub(jh) )
481	anaf(jh) = ftide ( ntide_sub(jh) )
482	END DO
483
484	fjulday_startharm=fjulday !Set this at very start and store
485
486	IF (lwp) THEN
487	WRITE(numout,*) '--------------------------'
488	WRITE(numout,*) ' - Output anaf for check'
489	WRITE(numout,*) 'ANA F', anaf
490	WRITE(numout,*) 'ANA U', anau
491	WRITE(numout,*) 'ANA V', anav
492	WRITE(numout,*) fjulday_startharm
493	WRITE(numout,*) '--------------------------'
494	ENDIF
495
496	ENDIF
497
498	ELSE
499
500	IF (lwp) WRITE(numout,*) "No variable setup for harmonic analysis"
501
502	ENDIF
503
504	END SUBROUTINE harm_ana_init
505	!
506	SUBROUTINE gelim (a,b,x,n)
507	!!----------------------------------------------------------------------
508	!! * ROUTINE harm_ana *
509	!!
510	!! ** Purpose : Guassian elimination
511	!!
512	!!
513	!! ** Action : - first action (share memory array/varible modified
514	!! in this routine
515	!! - second action .....
516	!! - .....
517	!!
518	!! References :
519	!! Give references if exist otherwise suppress these lines
520	!!
521	!! History :
522	implicit none
523	!
524	integer :: n
525	REAL(WP) :: b(nb_ana2+1), a(nb_ana2+1,nb_ana*2+1)
526	REAL(WP) :: x(nb_ana*2+1)
527	INTEGER :: row,col,prow,pivrow,rrow
528	REAL(WP) :: atemp
529	REAL(WP) :: pivot
530	REAL(WP) :: m
531
532	do row=1,n-1
533	pivrow=row
534	pivot=a(row,n-row+1)
535	do prow=row+1,n
536	if (abs(a(prow,n-row+1)).gt.abs(pivot) ) then
537	pivot=a(prow,n-row+1)
538	pivrow=prow
539	endif
540	enddo
541	! swap row and prow
542	if ( pivrow .ne. row ) then
543	atemp=b(pivrow)
544	b(pivrow)=b(row)
545	b(row)=atemp
546	do col=1,n
547	atemp=a(pivrow,col)
548	a(pivrow,col)=a(row,col)
549	a(row,col)=atemp
550	enddo
551	endif
552
553	do rrow=row+1,n
554	if (a(row,row).ne.0) then
555
556	m=-a(rrow,n-row+1)/a(row,n-row+1)
557	do col=1,n
558	a(rrow,col)=m*a(row,col)+a(rrow,col)
559	enddo
560	b(rrow)=m*b(row)+b(rrow)
561	endif
562	enddo
563	enddo
564	! back substitution now
565
566	x(1)=b(n)/a(n,1)
567	do row=n-1,1,-1
568	x(n-row+1)=b(row)
569	do col=1,(n-row)
570	x(n-row+1)=(x(n-row+1)-a(row,col)*x(col))
571	enddo
572
573	x(n-row+1)=(x(n-row+1)/a(row,(n-row)+1))
574	enddo
575
576	return
577	END SUBROUTINE gelim
578
579	SUBROUTINE harm_ana_out
580	!!----------------------------------------------------------------------
581	!! * ROUTINE harm_ana_init *
582	!!
583	!! ** Purpose : initialization of ....
584	!!
585	!! ** Method : blah blah blah ...
586	!!
587	!! ** input : Namlist namexa
588	!!
589	!! ** Action : ...
590	!!
591	!! history :
592	!! 9.0 ! 03-08 (Autor Names) Original code
593	!!----------------------------------------------------------------------
594	USE dianam ! build name of file (routine)
595
596	!! * local declarations
597	INTEGER :: ji, jj, jk, jgrid, jh ! dummy loop indices
598	! INTEGER :: nh_T
599	! INTEGER :: nid_harm
600	! CHARACTER (len=40) :: clhstnamt, clop1, clop2 ! temporary names
601	! CHARACTER (len=40) :: clhstnamu, clhstnamv ! temporary names
602	CHARACTER (len=40) :: suffix
603	! REAL(wp) :: zsto1, zsto2, zout, zmax, zjulian, zdt, zmdi ! temporary scalars
604
605	do jgrid=1,nvar_2d
606	do jh=1,nb_ana
607	h_out2D = 0.0
608	g_out2D = 0.0
609	do jj=1,nlcj
610	do ji=1,nlci
611	cca=g_cosamp2D(jh,ji,jj,jgrid)
612	ssa=g_sinamp2D(jh,ji,jj,jgrid)
613	h_out2D(ji,jj)=sqrt(cca2+ssa2)
614	IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
615	g_out2D(ji,jj)= 0.0_wp
616	ELSE
617	g_out2D(ji,jj)=(180.0/rpi)*atan2(ssa,cca)
618	ENDIF
619	IF (h_out2D(ji,jj).ne.0) THEN
620	h_out2D(ji,jj)=h_out2D(ji,jj)/anaf(jh)
621	ENDIF
622	IF (g_out2D(ji,jj).ne.0) THEN !Correct and take modulus
623	g_out2D(ji,jj) = g_out2D(ji,jj) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
624	if (g_out2D(ji,jj).gt.360.0) then
625	g_out2D(ji,jj)=g_out2D(ji,jj)-360.0
626	else if (g_out2D(ji,jj).lt.0.0) then
627	g_out2D(ji,jj)=g_out2D(ji,jj)+360.0
628	endif
629	ENDIF
630	enddo
631	enddo
632	!
633	! NETCDF OUTPUT
634	suffix = TRIM( m_varName2d( m_posi_2d(jgrid) ) )
635	CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix), h_out2D(:,:) )
636	CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix), g_out2D(:,:) )
637
638	enddo
639	enddo
640	!
641	! DO THE SAME FOR 3D VARIABLES
642	!
643	do jgrid=1,nvar_3d
644	do jh=1,nb_ana
645	h_out3D = 0.0
646	g_out3D = 0.0
647	DO jk=1,jpkm1
648	do jj=1,nlcj
649	do ji=1,nlci
650	cca=g_cosamp3D(jh,ji,jj,jk,jgrid)
651	ssa=g_sinamp3D(jh,ji,jj,jk,jgrid)
652	h_out3D(ji,jj,jk)=sqrt(cca2+ssa2)
653	IF (cca.eq.0.0 .and. ssa.eq.0.0) THEN
654	g_out3D(ji,jj,jk) = 0.0_wp
655	ELSE
656	g_out3D(ji,jj,jk) = (180.0/rpi)*atan2(ssa,cca)
657	ENDIF
658	IF (h_out3D(ji,jj,jk).ne.0) THEN
659	h_out3D(ji,jj,jk) = h_out3D(ji,jj,jk)/anaf(jh)
660	ENDIF
661	IF (g_out3D(ji,jj,jk).ne.0) THEN !Correct and take modulus
662	g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk) + MOD( (anau(jh)+anav(jh))/rad , 360.0)
663	if (g_out3D(ji,jj,jk).gt.360.0) then
664	g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)-360.0
665	else if (g_out3D(ji,jj,jk).lt.0.0) then
666	g_out3D(ji,jj,jk) = g_out3D(ji,jj,jk)+360.0
667	endif
668	ENDIF
669	enddo ! ji
670	enddo ! jj
671	ENDDO ! jk
672	!
673	! NETCDF OUTPUT
674	suffix = TRIM( m_varName3d( m_posi_3d(jgrid) ) )
675	IF(lwp) WRITE(numout,*) "harm_ana_out", suffix
676	CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'amp_'//TRIM(suffix), h_out3D(:,:,:) )
677	CALL iom_put( TRIM(Wave(ntide_all(jh))%cname_tide)//'pha_'//TRIM(suffix), g_out3D(:,:,:) )
678	enddo ! jh
679	enddo ! jgrid
680	!
681	END SUBROUTINE harm_ana_out
682	!
683	SUBROUTINE harm_rst_write(kt)
684	!!----------------------------------------------------------------------
685	!! * ROUTINE harm_ana_init *
686	!!
687	!! ** Purpose : To write out cummulated Tidal Harmomnic data to file for
688	!! restarting
689	!!
690	!! ** Method : restart files will be dated by default
691	!!
692	!! ** input :
693	!!
694	!! ** Action : ...
695	!!
696	!! history :
697	!! 0.0 ! 01-16 (Enda O'Dea) Original code
698	!! ASSUMES dated file for rose , can change later to be more generic
699	!!----------------------------------------------------------------------
700	INTEGER, INTENT(in) :: kt ! ocean time-step
701	!!
702	INTEGER :: jh, j2d, j3d
703	CHARACTER(LEN=20) :: clkt ! ocean time-step define as a character
704	CHARACTER(LEN=50) :: clname ! ocean output restart file name
705	CHARACTER(LEN=150) :: clpath ! full path to ocean output restart file
706	CHARACTER(LEN=250) :: clfinal ! full name
707
708	!restart file
709	DO j2d=1,nvar_2d
710	CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
711	DO jh=1,nb_ana
712	CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2 , :, :, j2d ) )
713	CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) )
714	ENDDO
715	ENDDO
716
717	DO j3d=1,nvar_3d
718	CALL iom_rstput( kt, nitrst, numrow, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
719	DO jh=1,nb_ana
720	CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2 , :, :, :, j3d ) )
721	CALL iom_rstput( kt, nitrst, numrow, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) )
722	ENDDO
723	ENDDO
724
725	IF(lwp) THEN
726	IF( kt > 999999999 ) THEN ; WRITE(clkt, * ) kt
727	ELSE ; WRITE(clkt, '(i8.8)') kt
728	ENDIF
729	clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
730	clpath = TRIM(cn_ocerst_outdir)
731	IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
732	IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for writing: ',TRIM(clpath)//clname
733
734	WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
735	OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
736	WRITE(66) cc
737	WRITE(66) anau
738	WRITE(66) anav
739	WRITE(66) anaf
740	WRITE(66) fjulday_startharm
741	CLOSE(66)
742	WRITE(numout,*) '----------------------------'
743	WRITE(numout,*) ' harm_rst_write: DONE '
744	WRITE(numout,*) cc
745	WRITE(numout,*) anaf
746	WRITE(numout,*) fjulday_startharm
747	WRITE(numout,*) '----------------------------'
748	ENDIF
749
750	END SUBROUTINE harm_rst_write
751
752	SUBROUTINE harm_rst_read
753	!!----------------------------------------------------------------------
754	!! * ROUTINE harm_ana_init *
755	!!
756	!! ** Purpose : To read in cummulated Tidal Harmomnic data to file for
757	!! restarting
758	!!
759	!! ** Method :
760	!!
761	!! ** input :
762	!!
763	!! ** Action : ...
764	!!
765	!! history :
766	!! 0.0 ! 01-16 (Enda O'Dea) Original code
767	!! ASSUMES dated file for rose , can change later to be more generic
768	!!----------------------------------------------------------------------
769	CHARACTER(LEN=20) :: clkt ! ocean time-step define as a character
770	CHARACTER(LEN=50) :: clname ! ocean output restart file name
771	CHARACTER(LEN=150) :: clpath ! full path to ocean output restart file
772	CHARACTER(LEN=250) :: clfinal ! full name
773	INTEGER :: jh, j2d, j3d
774
775	IF( nit000 > 999999999 ) THEN ; WRITE(clkt, * ) nit000-1
776	ELSE ; WRITE(clkt, '(i8.8)') nit000-1
777	ENDIF
778	clname = TRIM(cexper)//"_"//TRIM(ADJUSTL(clkt))//"_restart_harm_ana.bin"
779	clpath = TRIM(cn_ocerst_outdir)
780	IF( clpath(LEN_TRIM(clpath):) /= '/' ) clpath = TRIM(clpath) // '/'
781
782	IF (lwp) WRITE(numout,*) 'Open tidal harmonics restart file for reading: ',TRIM(clpath)//clname
783
784	DO j2d=1,nvar_2d
785	CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_2d(j2d) )), g_cumul_var2D( 1, :, :, j2d ) )
786	IF(lwp) WRITE(numout,*) "2D", j2d, m_posi_2d(j2d), m_varName2d( m_posi_2d(j2d) )
787	DO jh=1,nb_ana
788	CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_cos', g_cumul_var2D( jh*2 , :, :, j2d ) )
789	CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName2d( m_posi_2d(j2d) ))//'_sin', g_cumul_var2D( jh*2+1, :, :, j2d ) )
790	ENDDO
791	ENDDO
792
793	DO j3d=1,nvar_3d
794	CALL iom_get( numror,jpdom_autoglo, 'Mean_'//TRIM(m_varName2d( m_posi_3d(j3d) )), g_cumul_var3D( 1, :, :, :, j3d ) )
795	IF(lwp) WRITE(numout,*) "3D", j3d, m_posi_3d(j3d), m_varName3d( m_posi_3d(j3d) )
796
797	DO jh=1,nb_ana
798	CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_cos', g_cumul_var3D( jh*2 , :, :, :, j3d ) )
799	CALL iom_get( numror,jpdom_autoglo, TRIM(Wave(ntide_all(jh))%cname_tide)//"_"//TRIM(m_varName3d( m_posi_3d(j3d) ))//'_sin', g_cumul_var3D( jh*2+1, :, :, :, j3d ) )
800	ENDDO
801	ENDDO
802
803	WRITE(clfinal,'(a)') trim(clpath)//trim(clname)
804	OPEN( 66, file=TRIM(clfinal), form='unformatted', access="stream" )
805	READ(66) cc
806	READ(66) anau
807	READ(66) anav
808	READ(66) anaf
809	READ(66) fjulday_startharm
810	CLOSE(66)
811
812	IF(lwp) THEN
813	WRITE(numout,*) '----------------------------'
814	WRITE(numout,*) ' Checking anaf is correct'
815	WRITE(numout,*) cc
816	WRITE(numout,*) anaf
817	WRITE(numout,*) fjulday_startharm
818	WRITE(numout,*) '----------------------------'
819	ENDIF
820
821	END SUBROUTINE harm_rst_read
822
823	!!======================================================================
824	#else
825	!!---------------------------------------------------------------------------------
826	!! Dummy module NO harmonic Analysis
827	!!---------------------------------------------------------------------------------
828	LOGICAL, PUBLIC, PARAMETER :: lk_diaharm_fast = .FALSE. ! to be run or not
829
830	CONTAINS
831	SUBROUTINE harm_rst_write(kt) ! Dummy routine
832	END SUBROUTINE harm_rst_write
833	SUBROUTINE harm_rst_read ! Dummy routine
834	END SUBROUTINE harm_rst_read
835	SUBROUTINE harm_ana_out ! Dummy routine
836	END SUBROUTINE harm_ana_out
837	SUBROUTINE harm_ana_init
838	END SUBROUTINE harm_ana_init
839	SUBROUTINE harm_ana( kt )
840	!--- NB : end call not properly written
841	END SUBROUTINE harm_ana
842	! END SUBROUTINE harm_ana_init
843	!--- END NB
844	SUBROUTINE gelim (a,b,x,n)
845	!--- NB : end call not properly written
846	END SUBROUTINE gelim
847	! END SUBROUTINE gelim (a,b,x,n)
848	!--- END NB
849	#endif
850
851	END MODULE diaharm_fast

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

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

Context Navigation

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

Download in other formats: