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 : 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_mean, cdfid ! NETCDF file IDs.
38 integer :: id_var ! NETCDF variable ID.
39 integer :: ivtype ! 4=integer, 5=real, 6=d.p.
40 integer :: ndims ! Number of field dimensions.
41 integer :: natts ! Number of field attributes.
42 real :: member_inv ! 1 / ensemble size.
43
44 character(len=10) :: cv(1:max_num_vars) ! Default array of variable names.
45 integer :: dimids(1:10) ! Array of dimension IDs.
46 integer :: dims(1:max_num_dims) ! Array of dimensions.
47 integer :: istart(1:max_num_dims) ! Array of dimension starts.
48 integer :: iend(1:max_num_dims) ! Array of dimension ends.
49
50 real (kind=4), allocatable :: data_r(:,:,:) ! Data array.
51 real (kind=4), allocatable :: data_r_mean(:,:,:) ! Data array mean.
52
53 namelist / gen_be_ensmean_nl / filestub, num_members, nv, cv
54
55 stderr = 0
56 stdout = 6
57
58 !---------------------------------------------------------------------------------------------
59 write(6,'(/a)')' [1] Initialize information.'
60 !---------------------------------------------------------------------------------------------
61
62 filestub = 'test'
63 num_members = 56
64 nv = 1
65 cv = "U"
66
67 open(unit=unit, file='gen_be_ensmean_nl.nl', &
68 form='formatted', status='old', action='read')
69 read(unit, gen_be_ensmean_nl)
70 close(unit)
71
72 write(6,'(a,a)')' Filestub = ', trim(filestub)
73 write(6,'(a,i4)')' Number of ensemble members = ', num_members
74 write(6,'(a,i4)')' Number of variables to average = ', nv
75 write(6,'(50a)')' List of variables to average = ', cv(1:nv)
76
77 ! Open template ensemble mean with write access:
78 length = len_trim(filestub)
79 rcode = nf_open(filestub(1:length), NF_WRITE, cdfid_mean )
80 if ( rcode /= 0) then
81 write(UNIT=message(1),FMT='(A,A)') &
82 ' Error opening netcdf file ', filestub(1:length)
83 call da_error(__FILE__,__LINE__,message(1:1))
84 end if
85
86 !---------------------------------------------------------------------------------------------
87 write(6,'(/a)')' [4] Extract necessary fields from WRF ensemble forecasts.'
88 !---------------------------------------------------------------------------------------------
89
90 do v = 1, nv ! Loop over variables to average:
91 var = cv(v)
92 write(6,'(2a)')' Computing ensemble mean for variable ', var
93
94 do member = 1, num_members
95 write(UNIT=ce,FMT='(i3.3)')member
96
97 ! Open file:
98 input_file = trim(filestub)//'.e'//ce
99 length = len_trim(input_file)
100 rcode = nf_open( input_file(1:length), NF_NOWRITE, cdfid )
101
102 if ( member == 1 ) then
103 ! Get variable ID:
104 rcode = nf_inq_varid ( cdfid, var, id_var )
105
106 ! Check variable is in file:
107 if ( rcode /= 0 ) then
108 write(UNIT=message(1),FMT='(A,A)') &
109 var, ' variable is not in input file'
110 call da_error(__FILE__,__LINE__,message(1:1))
111 end if
112
113 ! Get dimension information for this field:
114 dimids = 0
115 dims = 0
116 rcode = nf_inq_var( cdfid, id_var, var, ivtype, ndims, dimids, natts )
117 do i = 1, ndims
118 rcode = nf_inq_dimlen( cdfid, dimids(i), dims(i) )
119 end do
120 istart = 1
121 iend = 1
122 do i = 1, ndims-1
123 iend(i) = dims(i)
124 end do
125
126 ! Allocate and initialize data:
127 if ( ivtype == 5 ) then
128 allocate( data_r(iend(1),iend(2),iend(3)))
129 allocate( data_r_mean(iend(1),iend(2),iend(3)))
130 data_r_mean = 0.0
131 else
132 write(UNIT=message(1),FMT='(A,A)') &
133 var, ' variable is not real type'
134 call da_error(__FILE__,__LINE__,message(1:1))
135 end if
136
137 end if
138
139 ! Calculate accumulating mean:
140 member_inv = 1.0 / real(member)
141 call ncvgt( cdfid, id_var, istart, iend, data_r, rcode)
142 data_r_mean = ( real( member-1 ) * data_r_mean + data_r ) * member_inv
143
144 rcode = nf_close( cdfid )
145
146 end do ! member
147
148 call ncvpt( cdfid_mean, id_var, istart, iend, data_r_mean, rcode)
149 deallocate( data_r )
150 deallocate( data_r_mean )
151
152 end do ! variable
153
154 rcode = nf_close( cdfid_mean )
155
156 end program gen_be_ensmean
157