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