gen_be_ensmean.f90
References to this file elsewhere.
1 program gen_be_ensmean
2 !
3 !----------------------------------------------------------------------
4 ! Purpose: Calculate NETCDF wrfinput format ensemble mean file from
5 ! input 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 : trace_use, stdout, stderr,filename_len
17 use da_tracing, only : da_trace_init, da_trace_entry, da_trace_exit, &
18 da_trace_report
19 use da_reporting, only : da_error,message
20
21 implicit none
22
23 #include <netcdf.inc>
24
25 integer, parameter :: max_num_dims = 4 ! Maximum number of dimensions.
26 integer, parameter :: max_num_vars = 50 ! Maximum number of variables.
27 integer, parameter :: unit = 100 ! Unit number.
28
29 character (len=filename_len) :: filestub ! General filename stub.
30 character (len=filename_len) :: input_file ! Input file.
31 character (len=10) :: var ! Variable to search for.
32 character (len=3) :: ce ! Member index -> character.
33
34 integer :: num_members ! Ensemble size.
35 integer :: nv ! Number of variables.
36 integer :: v, member, i ! Loop counters.
37 integer :: length ! Filename length.
38 integer :: rcode ! NETCDF return code.
39 integer :: cdfid_mean, cdfid ! 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
55 namelist / gen_be_ensmean_nl / filestub, num_members, nv, cv
56
57 stderr = 0
58 stdout = 6
59
60 !---------------------------------------------------------------------------------------------
61 write(6,'(/a)')' [1] Initialize information.'
62 !---------------------------------------------------------------------------------------------
63
64 if (trace_use) call da_trace_init
65 if (trace_use) call da_trace_entry("gen_be_ensmean")
66
67 filestub = 'test'
68 num_members = 56
69 nv = 1
70 cv = "U"
71
72 open(unit=unit, file='gen_be_ensmean_nl.nl', &
73 form='formatted', status='old', action='read')
74 read(unit, gen_be_ensmean_nl)
75 close(unit)
76
77 write(6,'(a,a)')' Filestub = ', trim(filestub)
78 write(6,'(a,i4)')' Number of ensemble members = ', num_members
79 write(6,'(a,i4)')' Number of variables to average = ', nv
80 write(6,'(50a)')' List of variables to average = ', cv(1:nv)
81
82 ! Open template ensemble mean with write access:
83 length = len_trim(filestub)
84 rcode = nf_open(filestub(1:length), NF_WRITE, cdfid_mean )
85 if ( rcode /= 0) then
86 write(UNIT=message(1),FMT='(A,A)') &
87 ' Error opening netcdf file ', filestub(1:length)
88 call da_error(__FILE__,__LINE__,message(1:1))
89 end if
90
91 !---------------------------------------------------------------------------------------------
92 write(6,'(/a)')' [4] Extract necessary fields from WRF ensemble forecasts.'
93 !---------------------------------------------------------------------------------------------
94
95 do v = 1, nv ! Loop over variables to average:
96 var = cv(v)
97 write(6,'(2a)')' Computing ensemble mean for variable ', var
98
99 do member = 1, num_members
100 write(UNIT=ce,FMT='(i3.3)')member
101
102 ! Open file:
103 input_file = trim(filestub)//'.e'//ce
104 length = len_trim(input_file)
105 rcode = nf_open( input_file(1:length), NF_NOWRITE, cdfid )
106
107 if ( member == 1 ) then
108 ! Get variable ID:
109 rcode = nf_inq_varid ( cdfid, var, id_var )
110
111 ! Check variable is in file:
112 if ( rcode /= 0 ) then
113 write(UNIT=message(1),FMT='(A,A)') &
114 var, ' variable is not in input file'
115 call da_error(__FILE__,__LINE__,message(1:1))
116 end if
117
118 ! Get dimension information for this field:
119 dimids = 0
120 dims = 0
121 rcode = nf_inq_var( cdfid, id_var, var, ivtype, ndims, dimids, natts )
122 do i = 1, ndims
123 rcode = nf_inq_dimlen( cdfid, dimids(i), dims(i) )
124 end do
125 istart = 1
126 iend = 1
127 do i = 1, ndims-1
128 iend(i) = dims(i)
129 end do
130
131 ! Allocate and initialize data:
132 if ( ivtype == 5 ) then
133 allocate( data_r(iend(1),iend(2),iend(3)))
134 allocate( data_r_mean(iend(1),iend(2),iend(3)))
135 data_r_mean = 0.0
136 else
137 write(UNIT=message(1),FMT='(A,A)') &
138 var, ' variable is not real type'
139 call da_error(__FILE__,__LINE__,message(1:1))
140 end if
141
142 end if
143
144 ! Calculate accumulating mean:
145 member_inv = 1.0 / real(member)
146 call ncvgt( cdfid, id_var, istart, iend, data_r, rcode)
147 data_r_mean = ( real( member-1 ) * data_r_mean + data_r ) * member_inv
148
149 rcode = nf_close( cdfid )
150
151 end do ! member
152
153 call ncvpt( cdfid_mean, id_var, istart, iend, data_r_mean, rcode)
154 deallocate( data_r )
155 deallocate( data_r_mean )
156
157 end do ! variable
158
159 rcode = nf_close( cdfid_mean )
160
161 if (trace_use) call da_trace_exit("gen_be_ensmean")
162 if (trace_use) call da_trace_report
163
164 end program gen_be_ensmean
165