gen_be_diags_read.f90
References to this file elsewhere.
1 program gen_be_diags_read
2
3 use da_control
4 use da_gen_be
5
6 implicit none
7
8 character*10 :: variable ! Variable name
9 character*8 :: uh_method ! Uh_method (power, scale)
10 character(len=filename_len) :: filename ! Input filename.
11 integer :: outunit ! Output unit for diagnostics.
12 integer :: ni, nj, nk, nk_3d ! Dimensions read in.
13 integer :: bin_type ! Type of bin to average over. !!!DALE ADD.
14 integer :: num_bins ! Number of 3D bins.
15 integer :: num_bins2d ! Number of 2D bins.
16 integer :: k ! Loop counter.
17 integer :: kdum ! Dummy vertical index.
18 integer :: max_wavenumber ! Smallest scale required (ni/2 - 1).
19 real :: lat_min, lat_max ! Used if bin_type = 2 (degrees).
20 real :: binwidth_lat ! Used if bin_type = 2 (degrees). !!!DALE ADD..
21 real :: binwidth_hgt ! Used if bin_type = 2 (m). !!!DALE ADD..
22 real :: hgt_min, hgt_max ! Used if bin_type = 2 (m).
23 real :: scale_length_ps_u ! Scale length for scalar ps_u.
24
25 integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point.
26 integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point.
27 real, allocatable :: regcoeff1(:) ! psi/chi regression cooefficient.
28 real, allocatable :: regcoeff2(:,:) ! psi/ps regression cooefficient.
29 real, allocatable :: regcoeff3(:,:,:) ! psi/T regression cooefficient.
30 real, allocatable :: e_vec(:,:) ! Domain-averaged eigenvectors.
31 real, allocatable :: e_val(:) ! Domain-averaged eigenvalues.
32 real, allocatable :: e_vec_loc(:,:,:) ! Latitudinally varying eigenvectors.
33 real, allocatable :: e_val_loc(:,:) ! Latitudinally varying eigenvalues.
34 real, allocatable :: total_power(:) ! Total Power spectrum.
35 real, allocatable :: scale_length(:) ! Scale length for regional application.
36
37 namelist / gen_be_diags_nl / uh_method
38
39 integer :: iunit,namelist_unit
40
41 ! Hardwire because of complicated incrementing of unit numbers writing a file
42 ! each time
43 integer, parameter :: ounit = 71
44
45 stderr = 0
46 stdout = 6
47
48 call da_get_unit(iunit)
49 call da_get_unit(namelist_unit)
50
51 open(unit=namelist_unit, file='gen_be_diags_nl.nl', &
52 form='formatted', status='old', action='read')
53 read(namelist_unit, gen_be_diags_nl)
54 close(namelist_unit)
55
56 filename = 'gen_be.dat'
57 print '("*** Unit=",i3,3X,"filename=",a40)',iunit, filename
58 open (iunit, file = filename, form='unformatted')
59
60 !----------------------------------------------------------------------------
61 ! [1] Gather regression coefficients.
62 !----------------------------------------------------------------------------
63
64 ! Read the dimensions:
65 read(iunit)ni, nj, nk
66 nk_3d = nk
67
68 allocate( bin(1:ni,1:nj,1:nk) )
69 allocate( bin2d(1:ni,1:nj) )
70
71 ! Read bin info:
72
73 read(iunit)bin_type
74 read(iunit)lat_min, lat_max, binwidth_lat
75 read(iunit)hgt_min, hgt_max, binwidth_hgt
76 read(iunit)num_bins, num_bins2d
77 read(iunit)bin(1:ni,1:nj,1:nk)
78 read(iunit)bin2d(1:ni,1:nj)
79
80 ! Read the regression coefficients:
81 allocate( regcoeff1(1:num_bins) )
82 allocate( regcoeff2(1:nk,1:num_bins2d) )
83 allocate( regcoeff3(1:nk,1:nk,1:num_bins2d) )
84
85 read(iunit)regcoeff1
86 read(iunit)regcoeff2
87 read(iunit)regcoeff3
88
89 outunit = ounit + 100
90 call da_print_be_stats_p( outunit, ni, nj, nk, num_bins, num_bins2d, &
91 bin, bin2d, regcoeff1, regcoeff2, regcoeff3 )
92
93 !----------------------------------------------------------------------------
94 ! [2] Gather vertical error eigenvectors, eigenvalues.
95 !----------------------------------------------------------------------------
96
97 read(iunit)variable
98 read(iunit)nk, num_bins2d
99
100 allocate( e_vec(1:nk,1:nk) )
101 allocate( e_val(1:nk) )
102 allocate( e_vec_loc(1:nk,1:nk,1:num_bins2d) )
103 allocate( e_val_loc(1:nk,1:num_bins2d) )
104
105 read(iunit)e_vec
106 read(iunit)e_val
107 read(iunit)e_vec_loc
108 read(iunit)e_val_loc
109 call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
110 e_vec, e_val, e_vec_loc, e_val_loc )
111 read(iunit)variable
112 read(iunit)nk, num_bins2d
113 read(iunit)e_vec
114 read(iunit)e_val
115 read(iunit)e_vec_loc
116 read(iunit)e_val_loc
117 call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
118 e_vec, e_val, e_vec_loc, e_val_loc )
119
120 read(iunit)variable
121 read(iunit)nk, num_bins2d
122 read(iunit)e_vec
123 read(iunit)e_val
124 read(iunit)e_vec_loc
125 read(iunit)e_val_loc
126 call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
127 e_vec, e_val, e_vec_loc, e_val_loc )
128
129 read(iunit)variable
130 read(iunit)nk, num_bins2d
131 read(iunit)e_vec
132 read(iunit)e_val
133 read(iunit)e_vec_loc
134 read(iunit)e_val_loc
135 call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
136 e_vec, e_val, e_vec_loc, e_val_loc )
137
138 deallocate( e_vec )
139 deallocate( e_val )
140 deallocate( e_vec_loc )
141 deallocate( e_val_loc )
142
143 if (uh_method /= 'power') then
144 ! 2d fields: ps_u, ps:
145
146 read(iunit)variable
147 read(iunit)nk, num_bins2d
148
149 allocate( e_vec(1:nk,1:nk) )
150 allocate( e_val(1:nk) )
151 allocate( e_vec_loc(1:nk,1:nk,1:num_bins2d) )
152 allocate( e_val_loc(1:nk,1:num_bins2d) )
153
154 read(iunit)e_vec
155 read(iunit)e_val
156 read(iunit)e_vec_loc
157 read(iunit)e_val_loc
158 call da_print_be_stats_v( outunit, variable, nk, num_bins2d, &
159 e_vec, e_val, e_vec_loc, e_val_loc )
160 deallocate( e_vec )
161 deallocate( e_val )
162 deallocate( e_vec_loc )
163 deallocate( e_val_loc )
164 end if
165 ! To assign the dimension nk for 3d fields:
166 nk = nk_3d
167
168 !----------------------------------------------------------------------------
169 if (uh_method == 'power') then
170 ! write(6,'(/a)') '[3] Gather horizontal error power spectra.'
171 !----------------------------------------------------------------------------
172
173 do k = 1, nk
174 read(iunit)variable
175 read(iunit)max_wavenumber, kdum
176 if ( k == 1 ) allocate( total_power(0:max_wavenumber) )
177 read(iunit)total_power(:)
178 call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
179 end do
180
181 do k = 1, nk
182 read(iunit)variable
183 read(iunit)max_wavenumber, kdum
184 read(iunit)total_power(:)
185 call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
186 end do
187
188 do k = 1, nk
189 read(iunit)variable
190 read(iunit)max_wavenumber, kdum
191 read(iunit)total_power(:)
192 call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
193 end do
194
195 do k = 1, nk
196 read(iunit)variable
197 read(iunit)max_wavenumber, kdum
198 read(iunit)total_power(:)
199 call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
200 end do
201
202 read(iunit)variable
203 read(iunit)max_wavenumber, kdum
204 read(iunit)total_power(:)
205 call da_print_be_stats_h_global( outunit, variable, k, max_wavenumber, total_power )
206
207 else if (uh_method == 'scale ') then
208
209 allocate (scale_length(1:nk))
210
211 ! psi:
212 read(iunit) variable
213 read(iunit) scale_length
214 call da_print_be_stats_h_regional( outunit, variable, nk, scale_length )
215
216 ! chi_u:
217 read(iunit) variable
218 read(iunit) scale_length
219 call da_print_be_stats_h_regional( outunit, variable, nk, scale_length )
220
221 ! t_u:
222 read(iunit) variable
223 read(iunit) scale_length
224 call da_print_be_stats_h_regional( outunit, variable, nk, scale_length )
225
226 ! rh:
227 read(iunit) variable
228 read(iunit) scale_length
229 call da_print_be_stats_h_regional( outunit, variable, nk, scale_length )
230
231 ! ps_u:
232 read(iunit) variable
233 read(iunit) scale_length_ps_u
234 write(6,'(3a,i5)')' Scale length for variable ', trim(variable), ' in unit ', outunit
235 write(outunit,'(a,i4,1pe15.5)')trim(variable), 1, scale_length_ps_u
236 outunit = outunit + 1
237 write(6,*)
238
239 deallocate (scale_length)
240
241 endif
242
243 close(iunit)
244
245 end program gen_be_diags_read