gen_be_cov3d.f90
References to this file elsewhere.
1 program gen_be_cov3d
2
3 use da_control, only : stderr, stdout, filename_len
4 use da_tools_serial, only : da_get_unit,da_advance_cymdh
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) :: dat_dir ! Input data directory.
13 character(len=filename_len) :: filename ! Input filename.
14 character*3 :: ce ! Member index -> character.
15 integer :: ni, nj, nk ! Grid dimensions.
16 integer :: i, j, k, 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 :: field1(:,:,:) ! Field 1.
32 real, allocatable :: field2(:,:,:) ! Field 2.
33 integer, allocatable:: bin(:,:,:) ! Bin assigned to each 3D point.
34 integer, allocatable:: bin2d(:,:) ! Bin assigned to each 2D point.
35 integer, allocatable:: bin_pts(:) ! Number of points in bin (3D fields).
36 real, allocatable :: covar(:) ! Covariance between input fields.
37 real, allocatable :: var(:) ! Autocovariance of field.
38
39 namelist / gen_be_cov3d_nl / start_date, end_date, interval, &
40 ne, variable1, variable2, dat_dir
41
42 integer :: ounit,iunit,namelist_unit
43
44 stderr = 0
45 stdout = 6
46
47 !---------------------------------------------------------------------------------------------
48 write(6,'(a)')' [1] Initialize namelist variables and other scalars.'
49 !---------------------------------------------------------------------------------------------
50
51 call da_get_unit(ounit)
52 call da_get_unit(iunit)
53 call da_get_unit(namelist_unit)
54
55
56 start_date = '2004030312'
57 end_date = '2004033112'
58 interval = 24
59 ne = 1
60 variable1 = 'psi'
61 variable2 = 'chi'
62 dat_dir = '/mmmtmp1/dmbarker'
63
64 open(unit=namelist_unit, file='gen_be_cov3d_nl.nl', &
65 form='formatted', status='old', action='read')
66 read(namelist_unit, gen_be_cov3d_nl)
67 close(namelist_unit)
68
69 read(start_date(1:10), fmt='(i10)')sdate
70 read(end_date(1:10), fmt='(i10)')edate
71 write(6,'(4a)')' Computing covariance for fields ', variable1 , ' and ', variable2
72 write(6,'(4a)') ' Time period is ', start_date, ' to ', end_date
73 write(6,'(a,i8,a)')' Interval between dates = ', interval, 'hours.'
74 write(6,'(a,i8)')' Number of ensemble members at each time = ', ne
75
76 date = start_date
77 cdate = sdate
78
79 !---------------------------------------------------------------------------------------------
80 write(6,'(2a)')' [2] Read fields, and calculate correlations'
81 !---------------------------------------------------------------------------------------------
82
83 first_time = .true.
84
85 do while ( cdate <= edate )
86 write(6,'(a,a)')' Processing data for date ', date
87
88 do member = 1, ne
89
90 write(ce,'(i3.3)')member
91
92 ! Read first field:
93 filename = trim(variable1)//'/'//date(1:10)
94 filename = trim(filename)//'.'//trim(variable1)//'.e'//ce
95 open (iunit, file = filename, form='unformatted')
96 read(iunit)ni, nj, nk
97
98 if ( first_time ) then
99 write(6,'(a,3i8)')' i, j, k dimensions are ', ni, nj, nk
100 allocate( bin(1:ni,1:nj,1:nk) )
101 allocate( bin2d(1:ni,1:nj) )
102 allocate( field1(1:ni,1:nj,1:nk) )
103 allocate( field2(1:ni,1:nj,1:nk) )
104
105 ! Read bin info:
106 filename = 'bin.data'
107 open (iunit+1, file = filename, form='unformatted')
108 read(iunit+1)bin_type
109 read(iunit+1)lat_min, lat_max, binwidth_lat
110 read(iunit+1)hgt_min, hgt_max, binwidth_hgt
111 read(iunit+1)num_bins, num_bins2d
112 read(iunit+1)bin(1:ni,1:nj,1:nk)
113 read(iunit+1)bin2d(1:ni,1:nj)
114 close(iunit+1)
115
116 allocate( bin_pts(1:num_bins) )
117 allocate( covar(1:num_bins) )
118 allocate( var(1:num_bins) )
119 bin_pts(:) = 0
120 covar(:) = 0.0
121 var(:) = 0.0
122 first_time = .false.
123 end if
124
125 read(iunit)field1
126 close(iunit)
127
128 ! Read second field:
129 filename = trim(variable2)//'/'//date(1:10)
130 filename = trim(filename)//'.'//trim(variable2)//'.e'//ce
131 open (iunit, file = filename, form='unformatted')
132 read(iunit)ni, nj, nk
133 read(iunit)field2
134 close(iunit)
135
136 ! Calculate covariances:
137
138 do k = 1, nk
139 do j = 1, nj
140 do i = 1, ni
141 b = bin(i,j,k)
142 bin_pts(b) = bin_pts(b) + 1
143 coeffa = 1.0 / real(bin_pts(b))
144 coeffb = real(bin_pts(b)-1) * coeffa
145 covar(b) = coeffb * covar(b) + coeffa * field1(i,j,k) * field2(i,j,k)
146 var(b) = coeffb * var(b) + coeffa * field2(i,j,k) * field2(i,j,k)
147 end do
148 end do
149 end do
150 end do ! End loop over ensemble members.
151
152 ! Calculate next date:
153 call da_advance_cymdh( date, interval, new_date )
154 date = new_date
155 read(date(1:10), fmt='(i10)')cdate
156 end do ! End loop over times.
157
158 filename = trim(variable1)//'.'//trim(variable2)//'.dat'
159 open (ounit, file = filename, status='unknown')
160
161 do k = 1, nk
162 do j = 1, nj
163 do i = 1, ni
164 b = bin(i,j,k) ! Take value from center of i dimension.
165 if ( var(b) /= 0.0 ) then
166 write(ounit,'(f16.8)')covar(b) / var(b)
167 else
168 write(ounit,'(f16.8)')0.0
169 end if
170 end do
171 end do
172 end do
173
174 end program gen_be_cov3d
175