gen_be_cov2d.f90
References to this file elsewhere.
1 program gen_be_cov2d
2
3 use da_control
4 use da_gen_be
5
6 implicit none
7
8 character*10 :: start_date, end_date ! Starting and ending dates.
9 character*10 :: date, new_date ! Current date (ccyymmddhh).
10 character*10 :: variable1 ! Variable name
11 character*10 :: variable2 ! Variable name
12 character(len=filename_len) :: filename ! Input filename.
13 character*3 :: ce ! Member index -> character.
14 integer :: ni, nj, nk, nkdum ! Grid dimensions.
15 integer :: i, j, member ! Loop counters.
16 integer :: b ! Bin marker.
17 integer :: sdate, cdate, edate ! Starting, current ending dates.
18 integer :: interval ! Period between dates (hours).
19 integer :: ne ! Number of ensemble members.
20 integer :: bin_type ! Type of bin to average over.
21 integer :: num_bins ! Number of bins (3D fields).
22 integer :: num_bins2d ! Number of bins (3D fields).
23 real :: lat_min, lat_max ! Used if bin_type = 2 (degrees).
24 real :: binwidth_lat ! Used if bin_type = 2 (degrees).
25 real :: hgt_min, hgt_max ! Used if bin_type = 2 (m).
26 real :: binwidth_hgt ! Used if bin_type = 2 (m).
27 real :: coeffa, coeffb ! Accumulating mean coefficients.
28 logical :: first_time ! True if first file.
29
30 real, allocatable :: latitude(:,:) ! Latitude (degrees, from south).
31 real, allocatable :: height(:,:,:) ! Height field.
32 real, allocatable :: field1(:,:) ! Field 1.
33 real, allocatable :: field2(:,:) ! Field 2.
34 integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point.
35 integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point.
36 integer, allocatable:: bin_pts2d(:) ! Number of points in bin (2D fields).
37 real, allocatable :: covar(:) ! Covariance between input fields.
38 real, allocatable :: var(:) ! Autocovariance of field.
39
40 namelist / gen_be_cov2d_nl / start_date, end_date, interval, &
41 ne, bin_type, &
42 lat_min, lat_max, binwidth_lat, &
43 hgt_min, hgt_max, binwidth_hgt, &
44 variable1, variable2
45
46 integer :: ounit,iunit,namelist_unit
47
48 stderr = 0
49 stdout = 6
50
51
52 !---------------------------------------------------------------------------------------------
53 write(6,'(a)')' [1] Initialize namelist variables and other scalars.'
54 !---------------------------------------------------------------------------------------------
55
56
57 call da_get_unit(ounit)
58 call da_get_unit(iunit)
59 call da_get_unit(namelist_unit)
60
61
62 start_date = '2004030312'
63 end_date = '2004033112'
64 interval = 24
65 ne = 1
66 bin_type = 1
67 lat_min = -90.0
68 lat_max = 90.0
69 binwidth_lat = 10.0
70 hgt_min = 0.0
71 hgt_max = 20000.0
72 binwidth_hgt = 1000.0
73 variable1 = 'ps_u'
74 variable2 = 'ps'
75
76 open(unit=namelist_unit, file='gen_be_cov2d_nl.nl', &
77 form='formatted', status='old', action='read')
78 read(namelist_unit, gen_be_cov2d_nl)
79 close(namelist_unit)
80
81 read(start_date(1:10), fmt='(i10)')sdate
82 read(end_date(1:10), fmt='(i10)')edate
83 write(6,'(4a)')' Computing covariance for fields ', variable1 , ' and ', variable2
84 write(6,'(4a)') ' Time period is ', start_date, ' to ', end_date
85 write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
86 write(6,'(a,i8)')' Number of ensemble members at each time = ', ne
87
88 date = start_date
89 cdate = sdate
90
91 !---------------------------------------------------------------------------------------------
92 write(6,'(2a)')' [2] Read fields, and calculate correlations'
93 !---------------------------------------------------------------------------------------------
94
95 first_time = .true.
96
97 do while ( cdate <= edate )
98 write(6,'(a,a)')' Processing data for date ', date
99
100 do member = 1, ne
101
102 write(ce,'(i3.3)')member
103
104 ! Read Full-fields:
105 filename = 'fullflds/'//date(1:10)
106 filename = trim(filename)//'.fullflds.e'//ce
107 open (iunit, file = filename, form='unformatted')
108
109 read(iunit)ni, nj, nk
110 if ( first_time ) then
111 write(6,'(a,3i8)')' i, j, k dimensions are ', ni, nj, nk
112 allocate( latitude(1:ni,1:nj) )
113 allocate( height(1:ni,1:nj,1:nk) )
114 allocate( bin(1:ni,1:nj,1:nk) )
115 allocate( bin2d(1:ni,1:nj) )
116 allocate( field1(1:ni,1:nj) )
117 allocate( field2(1:ni,1:nj) )
118 end if
119 read(iunit)latitude
120 read(iunit)height
121
122 ! Create and sort into bins:
123 call da_create_bins( ni, nj, nk, bin_type, num_bins, num_bins2d, bin, bin2d, &
124 lat_min, lat_max, binwidth_lat, &
125 hgt_min, hgt_max, binwidth_hgt, latitude, height )
126
127 close(iunit)
128
129 if ( first_time ) then
130 allocate( bin_pts2d(1:num_bins2d) )
131 allocate( covar(1:num_bins2d) )
132 allocate( var(1:num_bins2d) )
133 bin_pts2d(:) = 0
134 covar(:) = 0.0
135 var(:) = 0.0
136 first_time = .false.
137 end if
138
139 ! Read 2D first field:
140 filename = trim(variable1)//'/'//date(1:10)
141 filename = trim(filename)//'.'//trim(variable1)//'.e'//ce//'.01'
142 open (iunit, file = filename, form='unformatted')
143 read(iunit)ni, nj, nkdum
144 read(iunit)field1
145 close(iunit)
146
147 ! Read 2D second field:
148 filename = trim(variable2)//'/'//date(1:10)
149 filename = trim(filename)//'.'//trim(variable2)//'.e'//ce//'.01'
150 open (iunit, file = filename, form='unformatted')
151 read(iunit)ni, nj, nkdum
152 read(iunit)field2
153 close(iunit)
154
155 ! Calculate covariances:
156
157 do j = 1, nj
158 do i = 1, ni
159 b = bin2d(i,j)
160 bin_pts2d(b) = bin_pts2d(b) + 1
161 coeffa = 1.0 / real(bin_pts2d(b))
162 coeffb = real(bin_pts2d(b)-1) * coeffa
163 covar(b) = coeffb * covar(b) + coeffa * field1(i,j) * field2(i,j)
164 var(b) = coeffb * var(b) + coeffa * field2(i,j) * field2(i,j)
165 end do
166 end do
167 end do ! End loop over ensemble members.
168
169 ! Calculate next date:
170 call da_advance_cymdh( date, interval, new_date )
171 date = new_date
172 read(date(1:10), fmt='(i10)')cdate
173 end do ! End loop over times.
174
175 filename = trim(variable1)//'.'//trim(variable2)//'.dat'
176 open (ounit, file = filename, status='unknown')
177
178 do j = 1, nj
179 do i = 1, ni
180 b = bin2d(i,j)
181 write(ounit,'(f16.8)')covar(b) / var(b)
182 end do
183 end do
184
185 end program gen_be_cov2d