gen_be_ensmean.f90
References to this file elsewhere.
1 program gen_be_ensmean
2 !
3 !----------------------------------------------------------------------
4 ! Purpose: Calculate ensemble mean file from input WRF NETCDF input
5 ! ensemble members.
6 !
7 ! Owner: Dale Barker (NCAR/MMM) - WRF wrappper. Thanks to Cindy Bruyere.
8 ! Please acknowledge author/institute in work that uses this code.
9 !
10 !----------------------------------------------------------------------
11
12 #ifdef crayx1
13 #define iargc ipxfargc
14 #endif
15
16 use da_control, only : stdout, stderr,filename_len
17 use da_reporting, only : da_error,message
18
19 implicit none
20
21 #include <netcdf.inc>
22
23 integer, parameter :: max_num_dims = 4 ! Maximum number of dimensions.
24 integer, parameter :: max_num_vars = 50 ! Maximum number of variables.
25 integer, parameter :: unit = 100 ! Unit number.
26
27 character (len=filename_len) :: filestub ! General filename stub.
28 character (len=filename_len) :: input_file ! Input file.
29 character (len=10) :: var ! Variable to search for.
30 character (len=3) :: ce ! Member index -> character.
31
32 integer :: num_members ! Ensemble size.
33 integer :: nv ! Number of variables.
34 integer :: v, member, i ! Loop counters.
35 integer :: length ! Filename length.
36 integer :: rcode ! NETCDF return code.
37 integer :: cdfid ! NETCDF file IDs.
38 integer :: cdfid_mean ! NETCDF file IDs.
39 integer :: cdfid_vari ! NETCDF file IDs.
40 integer :: id_var ! NETCDF variable ID.
41 integer :: ivtype ! 4=integer, 5=real, 6=d.p.
42 integer :: ndims ! Number of field dimensions.
43 integer :: natts ! Number of field attributes.
44 real :: member_inv ! 1 / ensemble size.
45
46 character(len=10) :: cv(1:max_num_vars) ! Default array of variable names.
47 integer :: dimids(1:10) ! Array of dimension IDs.
48 integer :: dims(1:max_num_dims) ! Array of dimensions.
49 integer :: istart(1:max_num_dims) ! Array of dimension starts.
50 integer :: iend(1:max_num_dims) ! Array of dimension ends.
51
52 real (kind=4), allocatable :: data_r(:,:,:) ! Data array.
53 real (kind=4), allocatable :: data_r_mean(:,:,:) ! Data array mean.
54 real (kind=4), allocatable :: data_r_vari(:,:,:) ! Data array variance.
55
56 namelist / gen_be_ensmean_nl / filestub, num_members, nv, cv
57
58 stderr = 0
59 stdout = 6
60
61 !---------------------------------------------------------------------------------------------
62 write(6,'(/a)')' [1] Initialize information.'
63 !---------------------------------------------------------------------------------------------
64
65 filestub = 'test'
66 num_members = 56
67 nv = 1
68 cv = "U"
69
70 open(unit=unit, file='gen_be_ensmean_nl.nl', &
71 form='formatted', status='old', action='read')
72 read(unit, gen_be_ensmean_nl)
73 close(unit)
74
75 write(6,'(a,a)')' Filestub = ', trim(filestub)
76 write(6,'(a,i4)')' Number of ensemble members = ', num_members
77 write(6,'(a,i4)')' Number of variables to average = ', nv
78 write(6,'(50a)')' List of variables to average = ', cv(1:nv)
79
80 ! Open template ensemble mean with write access:
81 length = len_trim(filestub)
82 rcode = nf_open(filestub(1:length), NF_WRITE, cdfid_mean )
83 if ( rcode /= 0) then
84 write(UNIT=message(1),FMT='(A,A)') &
85 ' Error opening netcdf file ', filestub(1:length)
86 call da_error(__FILE__,__LINE__,message(1:1))
87 end if
88
89 ! Open template ensemble variance with write access:
90 input_file = trim(filestub)//'.vari'
91 length = len_trim(input_file)
92 rcode = nf_open(input_file(1:length), NF_WRITE, cdfid_vari )
93 if ( rcode /= 0) then
94 write(UNIT=message(1),FMT='(A,A)') &
95 ' Error opening netcdf file ', input_file(1:length)
96 call da_error(__FILE__,__LINE__,message(1:1))
97 end if
98
99 !---------------------------------------------------------------------------------------------
100 write(6,'(/a)')' [2] Extract necessary fields from WRF ensemble forecasts.'
101 !---------------------------------------------------------------------------------------------
102
103 do v = 1, nv ! Loop over variables to average:
104 var = cv(v)
105 write(6,'(2a)')' Computing ensemble mean for variable ', var
106
107 do member = 1, num_members
108 write(UNIT=ce,FMT='(i3.3)')member
109
110 ! Open file:
111 input_file = trim(filestub)//'.e'//ce
112 length = len_trim(input_file)
113 rcode = nf_open( input_file(1:length), NF_NOWRITE, cdfid )
114
115 if ( member == 1 ) then
116 ! Get variable ID:
117 rcode = nf_inq_varid ( cdfid, var, id_var )
118
119 ! Check variable is in file:
120 if ( rcode /= 0 ) then
121 write(UNIT=message(1),FMT='(A,A)') &
122 var, ' variable is not in input file'
123 call da_error(__FILE__,__LINE__,message(1:1))
124 end if
125
126 ! Get dimension information for this field:
127 dimids = 0
128 dims = 0
129 rcode = nf_inq_var( cdfid, id_var, var, ivtype, ndims, dimids, natts )
130 do i = 1, ndims
131 rcode = nf_inq_dimlen( cdfid, dimids(i), dims(i) )
132 end do
133 istart = 1
134 iend = 1
135 do i = 1, ndims-1
136 iend(i) = dims(i)
137 end do
138
139 ! Allocate and initialize data:
140 if ( ivtype == 5 ) then
141 allocate( data_r(iend(1),iend(2),iend(3)))
142 allocate( data_r_mean(iend(1),iend(2),iend(3)))
143 allocate( data_r_vari(iend(1),iend(2),iend(3)))
144 data_r_mean = 0.0
145 data_r_vari = 0.0
146 else
147 write(UNIT=message(1),FMT='(A,A)') &
148 var, ' variable is not real type'
149 call da_error(__FILE__,__LINE__,message(1:1))
150 end if
151
152 end if
153
154 ! Calculate accumulating mean and variance:
155 member_inv = 1.0 / real(member)
156 call ncvgt( cdfid, id_var, istart, iend, data_r, rcode)
157 data_r_mean = ( real( member-1 ) * data_r_mean + data_r ) * member_inv
158 data_r_vari = ( real( member-1 ) * data_r_vari + data_r * data_r ) * member_inv
159
160 rcode = nf_close( cdfid )
161
162 end do ! member
163
164 call ncvpt( cdfid_mean, id_var, istart, iend, data_r_mean, rcode)
165 call ncvpt( cdfid_vari, id_var, istart, iend, data_r_vari, rcode)
166 deallocate( data_r )
167 deallocate( data_r_mean )
168 deallocate( data_r_vari )
169
170 end do ! variable
171
172 rcode = nf_close( cdfid_mean )
173 rcode = nf_close( cdfid_vari )
174
175 end program gen_be_ensmean
176