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