Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

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

Download in other formats: