gen_spectra.f90
References to this file elsewhere.
1 program gen_spectra
2
3 use da_control
4 use da_spectral
5 use da_tracing
6
7 implicit none
8
9 integer, parameter :: iunit = 11
10
11 real, parameter :: tolerance = 1.0e-6
12
13 character(len=filename_len) :: filename
14 character*10 :: date, cvar, ctime
15 integer :: time, member ! Loop counters.
16 integer :: i, j, k, n, m ! Loop counters.
17 integer :: jts, jte ! start, end indices.
18 integer :: nt ! Number of times to process.
19 integer :: ne ! Number of ensemble members.
20 integer :: ni, nj, nk ! Dimensions read in.
21 integer :: max_wavenumber ! Smallest scale required (ni/2 - 1).
22 integer :: index_k, index_start, index_end ! Array indexer.
23 integer :: inc ! Jump between elements of vector in array.
24 integer :: lenr ! FFT array dimension (at least inc*(n-1)+1).
25 integer :: lensav ! wsave dimension (n+int(log(real(ni)))+4).
26 integer :: lenwrk ! Dimension of work array.
27 integer :: ier ! FFT error flag.
28 integer :: count ! Counter
29 integer :: index_m, index_j ! Array indices.
30 integer :: alp_size ! Ass. Leg. Pol. array size.
31 integer :: cv_size ! Complex control variable array size.
32
33 logical :: first_time ! True if first time through loop.
34 logical :: testing_spectral ! Test spectral transform.
35 real :: pi_over_180 ! pi / 180
36 real :: diff_rms ! RMS error measure.
37
38 real, allocatable :: ps_prime(:,:) ! Surface pressure perturbation.
39 real, allocatable :: t_prime(:,:,:) ! Temperature perturbation.
40 real, allocatable :: psi_prime(:,:,:) ! Streamfunction perturbation.
41 real, allocatable :: chi_prime(:,:,:) ! Velocity Potential perturbation.
42 real, allocatable :: rh_prime(:,:,:) ! Relative Humidity Perturbation.
43 real, allocatable :: height(:,:,:) ! Geopotential height.
44 real, allocatable :: latitude(:,:) ! Latitude (radians)
45 real, allocatable :: field(:,:) ! Gridpoint field to be decomposed.
46 real, allocatable :: field_out(:,:) ! Output test field.
47
48 real, allocatable :: lat(:) ! Latitude (radians).
49 real, allocatable :: sinlat(:) ! sine(latitude).
50 real, allocatable :: coslat(:) ! cosine(latitude).
51 real, allocatable :: int_wgts(:) ! Legendre integration weights.
52 real, allocatable :: lon(:) ! Longitude (radians).
53 real, allocatable :: sinlon(:) ! sine(longitude).
54 real, allocatable :: coslon(:) ! cosine(longitude).
55 real, allocatable :: wsave(:) ! Prime values for fast Fourier transform.
56 real, allocatable :: alp(:) ! Associated Legendre Polynomial.
57 real, allocatable :: power(:) ! Power spectrum (n).
58 complex, allocatable:: cv(:) ! Control variable vector.
59 complex, allocatable:: cv_out(:) ! Control variable vector copy for tests.
60
61 if (trace_use) call da_trace_init
62 if (trace_use) call da_trace_entry("gen_spectra")
63
64 pi = 4.0 * atan(1.0)
65 pi_over_180 = pi / 180.0
66 gaussian_lats = .false.
67 testing_spectral = .true.
68 nt = 1
69 ne = 1
70
71 !---------------------------------------------------------------------------------------------
72 ! gen_statistics2: Read standard fields, calculate regression coefficients, and create and
73 ! output "control variables".
74 !---------------------------------------------------------------------------------------------
75
76 ! Loop over times, ensemble members.
77
78 !---------------------------------------------------------------------------------------------
79 write(6,'(a)')' [1] Read in perturbation fields.'
80 !---------------------------------------------------------------------------------------------
81
82 open (iunit, file='/data2/hcshin/youn/DIFF63/GHMD_DIFF.2004030212_old', form='unformatted')
83
84 read(iunit)date, ni, nj, nk
85 write(6,'(a,a10)')' Date = ', date
86 write(6,'(a,3i8)')' i, j, k dimensions are ', ni, nj, nk
87
88 allocate( ps_prime(1:ni,1:nj) )
89 allocate( t_prime(1:ni,1:nj,1:nk) )
90 allocate( psi_prime(1:ni,1:nj,1:nk) )
91 allocate( chi_prime(1:ni,1:nj,1:nk) )
92 allocate( rh_prime(1:ni,1:nj,1:nk) )
93 allocate( height(1:ni,1:nj,1:nk) )
94 allocate( latitude(1:ni,1:nj) )
95
96 read(iunit)ps_prime
97 read(iunit)t_prime
98 read(iunit)psi_prime
99 read(iunit)chi_prime
100 read(iunit)rh_prime
101 read(iunit)height
102 read(iunit)latitude ! Not currently used.
103
104 filename = "psi"//"date"//"member"//"mode"
105 open (iunit, file=filename, form='unformatted')
106 write(iunit)ni, nj
107 write(iunit)gaussian_lats
108 write(iunit)psi_prime(1:ni,1:nj,1)
109 close(iunit)
110
111 ! Calculate regression coefficients.
112
113 ! Calculate control variables.
114
115 ! Output vp(i,j,k) = vp^-1 vs(i,j,k).
116
117 !---------------------------------------------------------------------------------------------
118 ! gen_statistics3: Calculate Vertical Eigenvectors, Eigenvalues for a given variable.
119 !---------------------------------------------------------------------------------------------
120
121 ! Loop-over/distribute control variables.
122
123 ! Loop over times, ensemble members.
124
125 ! Read control variables
126
127 ! Calculate/accumulate Bv(j,k1,k2).
128
129 ! End loop over times, ensemble members.
130
131 ! Calculate E, L (v,j,k,k).
132
133 ! Loop over times, ensemble members:
134 ! Output v(i,j,m) = vv^-1 vp(i,j,k)
135
136 !---------------------------------------------------------------------------------------------
137 ! gen_statistics4: Calculate Power spectra for given 2D field.
138 !---------------------------------------------------------------------------------------------
139
140 ! Loop over/distribute variables, vertical modes:
141
142 ! Loop over times, ensemble members:
143
144 cvar = 'psi'
145
146 first_time = .true.
147
148 do time = 1, nt
149 write(ctime,'(i4)')time
150
151 do member = 1, ne
152 ! write(cmember,'(i4)')member
153
154
155 ! Read 2D fields.
156 filename = trim(cvar)//"date"//"member"//"mode"
157 open (iunit, file=filename, form='unformatted')
158 read(iunit)ni, nj
159 read(iunit)gaussian_lats
160
161 if ( first_time ) then ! Initialize
162 allocate( field(1:ni,1:nj) )
163 read(iunit)field
164
165 !---------------------------------------------------------------------------------------------
166 write(6,'(a)')' Initialize spectral transforms.'
167 !---------------------------------------------------------------------------------------------
168
169 inc = 1
170 lenr = inc * (ni - 1 ) + 1
171 lensav = ni + int(log(real(ni))) + 4
172 lenwrk = ni
173
174 allocate( lat(1:nj) )
175 allocate( sinlat(1:nj) )
176 allocate( coslat(1:nj) )
177 allocate( int_wgts(1:nj) )
178 allocate( lon(1:ni) )
179 allocate( sinlon(1:ni) )
180 allocate( coslon(1:ni) )
181 allocate( wsave(1:lensav) )
182
183 max_wavenumber = ni / 2 - 1
184 allocate ( power(0:max_wavenumber) )
185 power(:) = 0.0
186
187 alp_size = ( nj + 1 ) * ( max_wavenumber + 1 ) * ( max_wavenumber + 2 ) / 4
188 allocate( alp( 1:alp_size) )
189
190 call da_initialize_h( ni, nj, max_wavenumber, lensav, alp_size, &
191 wsave, lon, sinlon, coslon, lat, sinlat, coslat, &
192 int_wgts, alp )
193
194 cv_size = ( max_wavenumber + 1 ) * ( max_wavenumber + 2 ) / 2
195 allocate( cv( 1:cv_size) )
196
197 ! Test horizontal transforms:
198 if ( testing_spectral ) then
199 call da_test_spectral( ni, nj, max_wavenumber, inc, lenr, lensav, lenwrk, &
200 alp_size, cv_size, alp, wsave, int_wgts, field )
201 end if
202 first_time = .false.
203 else ! Just read field
204 read(iunit)field
205 end if
206
207 !---------------------------------------------------------------------------------------------
208 write(6,'(a)')' Perform gridpoint to spectral decomposition.'
209 !---------------------------------------------------------------------------------------------
210 call da_vv_to_v_spectral( ni, nj, max_wavenumber, inc, lenr, lensav, lenwrk, &
211 alp_size, cv_size, alp, wsave, int_wgts, field, cv )
212
213 !---------------------------------------------------------------------------------------------
214 write(6,'(a)')' Calculate power spectra.'
215 !---------------------------------------------------------------------------------------------
216
217 call da_calc_power( max_wavenumber, cv_size, cv, power )
218 do n = 0, max_wavenumber
219 write(6,'(i4,1pe15.5)')n, power(n)
220 end do
221
222 end do ! End loop over ensemble members.
223 end do ! End loop over times.
224
225 if (trace_use) call da_trace_exit("gen_spectra")
226 if (trace_use) call da_trace_report
227
228 end program gen_spectra