gen_be_stage4_global.f90
References to this file elsewhere.
1 program gen_be_stage4_global
2
3 use da_control, only : stderr, stdout, filename_len, pi, gaussian_lats
4 use da_be_spectral, only : da_vv_to_v_spectral, da_initialize_h, &
5 da_calc_power
6 use da_tools_serial, only : da_get_unit,da_advance_cymdh
7
8 implicit none
9
10 character*10 :: start_date, end_date ! Starting and ending dates.
11 character*10 :: date, new_date ! Current date (ccyymmddhh).
12 character*10 :: variable ! Variable name
13 character(len=filename_len) :: filename ! Input filename.
14 character*2 :: ck ! Loop index -> character.
15 character*3 :: ce ! Member index -> character.
16 integer :: ni, nj, nk ! Dimensions read in.
17 integer :: num_states ! Number of data times.
18 integer :: sdate, cdate, edate ! Starting, current ending dates.
19 integer :: interval ! Period between dates (hours).
20 integer :: k, member, n ! Loop counters.
21 integer :: ne ! Number of ensemble members.
22 integer :: max_wavenumber ! Smallest scale required (ni/2 - 1).
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 :: alp_size ! Ass. Leg. Pol. array size.
28 integer :: r_cvsize ! Real control variable array size.
29
30 logical :: first_time ! True if first time through loop.
31 logical :: testing_spectral ! True if testing spectral transforms.
32 real :: pi_over_180 ! pi / 180
33 real :: coeffa, coeffb ! Accumulating mean coefficients.
34 real :: variance ! Variance (sum of power spectra) .
35
36 real, allocatable :: field(:,:) ! Gridpoint field to be decomposed.
37
38 real, allocatable :: lat(:) ! Latitude (radians).
39 real, allocatable :: sinlat(:) ! sine(latitude).
40 real, allocatable :: coslat(:) ! cosine(latitude).
41 real, allocatable :: int_wgts(:) ! Legendre integration weights.
42 real, allocatable :: lon(:) ! Longitude (radians).
43 real, allocatable :: sinlon(:) ! sine(longitude).
44 real, allocatable :: coslon(:) ! cosine(longitude).
45 real, allocatable :: wsave(:) ! Prime values for fast Fourier transform.
46 real, allocatable :: alp(:) ! Associated Legendre Polynomial.
47 real, allocatable :: power(:) ! Power spectrum (n).
48 real, allocatable :: total_power(:) ! Total Power spectrum (averaged over time/members).
49
50 real, allocatable :: rcv(:) ! Control variable vector.
51
52 namelist / gen_be_stage4_global_nl / start_date, end_date, interval, variable, gaussian_lats, &
53 testing_spectral, ne, k
54
55 integer :: ounit,iunit,namelist_unit
56
57 stderr = 0
58 stdout = 6
59
60 call da_get_unit(ounit)
61 call da_get_unit(iunit)
62 call da_get_unit(namelist_unit)
63
64 pi_over_180 = pi / 180.0
65
66 write(unit=6,fmt='(a/)') "[1] Read in 2D perturbation fields for variable , variable"
67
68 start_date = '2004030312'
69 end_date = '2004033112'
70 interval = 24
71 variable = 'psi'
72 gaussian_lats = .false.
73 ne = 1
74 k = 1
75
76 open(unit=namelist_unit, file='gen_be_stage4_global_nl.nl', &
77 form='formatted', status='old', action='read')
78 read(namelist_unit, gen_be_stage4_global_nl)
79 close(namelist_unit)
80
81 write(ck,'(i2)')k
82 if ( k < 10 ) ck = '0'//ck(2:2)
83
84 read(start_date(1:10), fmt='(i10)')sdate
85 read(end_date(1:10), fmt='(i10)')edate
86 write(unit=6,fmt='(4a)') 'Computing horizontal power spectra for period' , start_date, 'to' , end_date
87 write(unit=6,fmt='(a,i8,a)') 'Interval between dates =' , interval,' hours'
88 write(unit=6,fmt='(a,i8)') 'Number of ensemble members at each time =', ne
89 write(unit=6,fmt='(3a,i4)') 'Variable' , variable,' at level' , k
90
91 date = start_date
92 cdate = sdate
93 first_time = .true.
94 num_states = 1
95
96 do while ( cdate <= edate )
97 do member = 1, ne
98 write(ce,'(i3.3)')member
99
100 ! write(6,(5a,i4)) Calculate spectra for date , date, , variable , trim(variable), &
101 ! and member , member
102
103 ! Read in data for given variable/level/time/member:
104
105 filename = trim(variable)//'/'//date(1:10)//'.'//trim(variable)
106 filename = trim(filename)//'.e'//ce//'.'//ck
107 open (iunit, file = filename, form='unformatted')
108 read(iunit)ni, nj, nk ! nk not used.
109
110 if ( first_time ) then
111 write(6,'(a,3i8)')' i, j, k dimensions are ', ni, nj, nk
112 allocate( field(1:ni,1:nj) )
113 end if
114 read(iunit)field
115 close(iunit)
116
117 if ( first_time ) then
118
119 write(unit=6,fmt='(a)') "Initialize spectral transforms"
120
121 inc = 1
122 lenr = inc * (ni - 1 ) + 1
123 lensav = ni + int(log(real(ni))) + 4
124 lenwrk = ni
125
126 allocate( lat(1:nj) )
127 allocate( sinlat(1:nj) )
128 allocate( coslat(1:nj) )
129 allocate( int_wgts(1:nj) )
130 allocate( lon(1:ni) )
131 allocate( sinlon(1:ni) )
132 allocate( coslon(1:ni) )
133 allocate( wsave(1:lensav) )
134
135 max_wavenumber = ni / 2 - 1
136 allocate ( power(0:max_wavenumber) )
137 allocate ( total_power(0:max_wavenumber) )
138 total_power(:) = 0.0
139
140 alp_size = ( nj + 1 ) * ( max_wavenumber + 1 ) * ( max_wavenumber + 2 ) / 4
141 allocate( alp( 1:alp_size) )
142
143 call da_initialize_h( ni, nj, max_wavenumber, lensav, alp_size, &
144 wsave, lon, sinlon, coslon, lat, sinlat, coslat, &
145 int_wgts, alp )
146
147 r_cvsize = ( max_wavenumber + 1 ) * ( max_wavenumber + 2 )
148 allocate( rcv( 1:r_cvsize) )
149
150 ! Test horizontal transforms:
151 ! if ( testing_spectral ) then
152 ! call da_test_spectral( ni, nj, max_wavenumber, inc, lenr, lensav, lenwrk, &
153 ! alp_size, r_cvsize, alp, wsave, int_wgts, field )
154 ! end if
155 first_time = .false.
156 end if
157
158 write(unit=6,fmt='(a)') "Perform gridpoint to spectral decomposition"
159
160 call da_vv_to_v_spectral( ni, nj, max_wavenumber, inc, lenr, lensav, lenwrk, &
161 alp_size, r_cvsize, alp, wsave, int_wgts, rcv, field )
162
163 write(unit=6,fmt='(a)') "Calculate power spectra"
164
165 call da_calc_power( max_wavenumber, r_cvsize, rcv, power )
166
167 coeffa = 1.0 / real(num_states)
168 coeffb = real(num_states-1) * coeffa
169
170 do n = 0, max_wavenumber
171 total_power(n) = total_power(n) * coeffb + power(n) * coeffa
172 ! write(6,(2i4,2f18.6))num_states, n, power(n), total_power(n)
173 end do
174
175 num_states = num_states + 1
176 end do ! End loop over ensemble members.
177
178 ! Calculate next date:
179 call da_advance_cymdh( date, interval, new_date )
180 date = new_date
181 read(date(1:10), fmt='(i10)')cdate
182
183 end do ! End loop over times.
184
185 variance = sum( total_power(0:max_wavenumber) )
186 write(6,'(3a,i2,a,1pe15.5)')' Variable = ', trim(variable), ', Vertical Index = ', &
187 k, ', Variance = ', variance
188
189 filename = trim(variable)//'/'//trim(variable)
190 filename = trim(filename)//'.'//ck//'.spectrum'
191 open (ounit, file = filename, form='unformatted')
192 write(ounit)variable
193 write(ounit)max_wavenumber, k
194 write(ounit)total_power
195
196 end program gen_be_stage4_global