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