1 | PROGRAM fbgenerate |
---|
2 | !!--------------------------------------------------------------------- |
---|
3 | !! |
---|
4 | !! ** PROGRAM fbmatchup ** |
---|
5 | !! |
---|
6 | !! ** Purpose : Generate a feedback file of pseudo obs |
---|
7 | !! |
---|
8 | !! ** Method : Use of utilities from obs_fbm. |
---|
9 | !! |
---|
10 | !! ** Action : Read in data from a namelist file and generate a feedback |
---|
11 | !! file of pseudo observations |
---|
12 | !! |
---|
13 | !! Usage: For any parameter, if 3 values are given, these will be treated as |
---|
14 | !! bounds (start, stop, step), EXCEPT when only 3 values are expected. |
---|
15 | !! |
---|
16 | !! For spatial coord, date and time |
---|
17 | !! 3 values with nlat=3 will be treated as three separate values |
---|
18 | !! 3 values with nobs=3 will be treated as three separate values |
---|
19 | !! 3 values with nobs>3 will be treated as a start, end and step |
---|
20 | !! Giving a FillValue in one of namelist entries may have unexpected consequences |
---|
21 | !! e.g., discarding rest of list of values |
---|
22 | !! |
---|
23 | !! For time bounds, start(HHMM), end(HHMM), step(minutes) |
---|
24 | !! For date bounds, start(YYYYMMDD), end(YYYYMMDD), step(days) |
---|
25 | !! |
---|
26 | !! To split obs evenly in time, date and time should be given two values |
---|
27 | !! describing the start and end times. |
---|
28 | !! |
---|
29 | !! Can use logical flag to shuffle timestamps to remove position-time correlation. |
---|
30 | !! |
---|
31 | !! Can use logical flag to add uniformly sampled perturbation (up to defined limit) to positions |
---|
32 | !! |
---|
33 | !! For obs values, can list either single value for all (x,y,z) |
---|
34 | !! or can specify single profile for all (x,y) |
---|
35 | !! or can specify all values 1 profile after another. |
---|
36 | !! |
---|
37 | !! In namelist to set vals(nlev,nobs,nvar) |
---|
38 | !! set all obs (for variable 1) at all levels to one value |
---|
39 | !! vals(:,:,1)= 1.0 |
---|
40 | !! set all obs (for variable 2) at all levels to one value |
---|
41 | !! vals(:,:,2)= 35.0 |
---|
42 | !! set one profile (for variable 1), apply to all |
---|
43 | !! vals(:,:,1)= 1.0, 2.0, 3.0, 4.0 |
---|
44 | !! set all profiles profile (for variable 1) |
---|
45 | !! vals(:,1,1)= 1.0, 2.0, 3.0, 4.0 |
---|
46 | !! vals(:,2,1)= 1.0, 2.0, 3.0, 4.0 |
---|
47 | !! vals(:,nobs,1)= 1.0, 2.0, 3.0, 4.0 |
---|
48 | !! |
---|
49 | !! |
---|
50 | !! Limitations: Uses an allocatable array in a namelist. This is not meant |
---|
51 | !! to be allowed in F95, but it works when compiled with xlf90_r |
---|
52 | !! (and is allowed in Fortran2003). |
---|
53 | !! |
---|
54 | !! Forces same depth levels on all profiles. |
---|
55 | !! |
---|
56 | !! History : |
---|
57 | !! ! 2014 (J. Waters) Initial version |
---|
58 | !! ! 2014-02 (R. King) Adapted for profiles and bounds on vars |
---|
59 | !!---------------------------------------------------------------------- |
---|
60 | USE obs_fbm |
---|
61 | USE fbgenerate_coords |
---|
62 | USE test_fbgenerate |
---|
63 | IMPLICIT NONE |
---|
64 | TYPE(obfbdata) :: fbobsdata |
---|
65 | CHARACTER(len=256) :: filenameout |
---|
66 | #ifndef NOIARGCPROTO |
---|
67 | INTEGER,EXTERNAL :: iargc |
---|
68 | #endif |
---|
69 | INTEGER :: nobs, nvar, nlev, ntyp, nadd, next, lgrid, ierr, nlats, nlons |
---|
70 | INTEGER :: i, j |
---|
71 | INTEGER :: FillValue_int |
---|
72 | REAL(KIND=fbdp) :: FillValue_real |
---|
73 | INTEGER, ALLOCATABLE :: dates(:), times(:) |
---|
74 | REAL(KIND=fbdp), ALLOCATABLE :: julian_dates(:) |
---|
75 | REAL(KIND=fbdp) :: phys_spacing, lat_step, current_lat, perturb_limit |
---|
76 | REAL(KIND=fbdp), ALLOCATABLE :: lats(:), lons(:), deps(:), vals(:,:,:), model_vals(:,:,:) |
---|
77 | CHARACTER(LEN=ilenwmo) :: cdwmo, cdtyp |
---|
78 | CHARACTER(LEN=ilenwmo), ALLOCATABLE :: variable(:) |
---|
79 | CHARACTER(LEN=ilenlong), ALLOCATABLE :: variable_longname(:) |
---|
80 | CHARACTER(LEN=ilenunit), ALLOCATABLE :: variable_units(:) |
---|
81 | LOGICAL, PARAMETER :: test_prog = .FALSE. |
---|
82 | LOGICAL :: ln_listobs, ln_physpace, ln_gridobs, ln_shuffletimes, ln_perturb_posn |
---|
83 | REAL(KIND=fbdp), PARAMETER :: earth_radius = 6371.0_fbdp |
---|
84 | REAL(KIND=fbdp), PARAMETER :: pi = 3.14159_fbdp |
---|
85 | |
---|
86 | !!!! set up default values and read in namelist namfbgen.in !!!!!!!!!!!!!!!!!!!! |
---|
87 | NAMELIST/namfbgeneratesetup/ln_listobs,ln_physpace,ln_gridobs,ln_shuffletimes,ln_perturb_posn,perturb_limit,nobs,nlats,nlons,phys_spacing,nlev,nvar,nadd |
---|
88 | NAMELIST/namfbgenerate/ntyp,lats,lons,deps,times,dates,vals,model_vals,variable,variable_longname,variable_units,cdwmo,cdtyp |
---|
89 | |
---|
90 | IF (test_prog) THEN |
---|
91 | WRITE(*,*) "In test mode:" |
---|
92 | CALL tester |
---|
93 | ELSE |
---|
94 | |
---|
95 | |
---|
96 | IF (iargc()/=1) THEN |
---|
97 | WRITE(*,*)'Usage:' |
---|
98 | WRITE(*,*)'fbgenerate <output filename>' |
---|
99 | CALL abort |
---|
100 | ENDIF |
---|
101 | |
---|
102 | CALL getarg(1,filenameout) |
---|
103 | |
---|
104 | ! set up default values for namfbgeneratesetup |
---|
105 | ln_listobs=.FALSE. |
---|
106 | ln_physpace=.FALSE. |
---|
107 | ln_gridobs=.FALSE. |
---|
108 | ln_shuffletimes=.FALSE. |
---|
109 | ln_perturb_posn=.FALSE. |
---|
110 | nobs=0 |
---|
111 | nlats=0 |
---|
112 | nlons=0 |
---|
113 | phys_spacing =100.0_fbdp |
---|
114 | nlev=0 |
---|
115 | nvar=0 |
---|
116 | |
---|
117 | ! Read in number of obs and number of depths from namelist (before reading other variables) |
---|
118 | FillValue_real = 99999.0_fbdp |
---|
119 | FillValue_int = 99999 |
---|
120 | OPEN(10,file='namfbgen.in') |
---|
121 | READ(10,namfbgeneratesetup) |
---|
122 | CLOSE(10) |
---|
123 | |
---|
124 | ! Calculate number of observations |
---|
125 | IF ((ln_gridobs) .AND. (.NOT.ln_physpace) .AND. (.NOT.ln_listobs)) THEN |
---|
126 | IF( ( nlats .LE. 0 ) .OR. (nlons .LE. 0) )THEN |
---|
127 | WRITE(6,*)'ERROR: nlats and nlons must be greater than 0 if ln_gridobs set' |
---|
128 | ELSE |
---|
129 | nobs=nlats*nlons |
---|
130 | WRITE(*,*) "Constructing a grid of lat-lons" |
---|
131 | ENDIF |
---|
132 | |
---|
133 | ELSE IF ((.NOT.ln_gridobs) .AND. (ln_physpace) .AND. (.NOT.ln_listobs)) THEN |
---|
134 | IF ((phys_spacing .LE. 0) .OR.(phys_spacing .GE. earth_radius)) THEN |
---|
135 | WRITE(6,*)'ERROR: phys_spacing must be greater than 0 and less than 6371km.' |
---|
136 | ELSE |
---|
137 | nlats = INT(pi * earth_radius / phys_spacing) |
---|
138 | IF (MOD(nlats,2)==0) nlats=nlats-1 ! Ensure nlats is odd to have even number around equator. |
---|
139 | lat_step = 180.0_fbdp / (nlats-1) |
---|
140 | nobs=0 |
---|
141 | current_lat=0.0_fbdp |
---|
142 | ! Sum up number of lons at each latitude |
---|
143 | DO i=1,(nlats-1)/2 |
---|
144 | current_lat = current_lat + lat_step * (pi/180.0_fbdp) |
---|
145 | nobs = nobs + INT(2.0_fbdp * pi * earth_radius * cos(current_lat) / phys_spacing) |
---|
146 | END DO |
---|
147 | nobs = nobs * 2 + INT(2.0_fbdp * pi * earth_radius / phys_spacing) |
---|
148 | WRITE(*,*) "Using a physical seperation of ~", INT(phys_spacing),"km" |
---|
149 | END IF |
---|
150 | |
---|
151 | ELSE IF ((.NOT.ln_gridobs) .AND. (.NOT.ln_physpace) .AND. (ln_listobs)) THEN |
---|
152 | IF( nobs .LE. 0 ) THEN |
---|
153 | WRITE(6,*)'ERROR: nobs must be greater than 0' |
---|
154 | ENDIF |
---|
155 | WRITE(*,*) "Constructing a list of lat-lons" |
---|
156 | ELSE |
---|
157 | WRITE(*,*)'ERROR: one (and only one) of the logical flags must be true!' |
---|
158 | CALL abort |
---|
159 | ENDIF |
---|
160 | |
---|
161 | |
---|
162 | ALLOCATE(dates(nobs), & |
---|
163 | times(nobs), & |
---|
164 | julian_dates(nobs), & |
---|
165 | lats(nobs), & |
---|
166 | lons(nobs), & |
---|
167 | deps(nlev), & |
---|
168 | variable(nvar), & |
---|
169 | variable_longname(nvar), & |
---|
170 | variable_units(nvar),& |
---|
171 | vals(nlev,nobs,nvar),& |
---|
172 | model_vals(nlev,nobs,nvar),& |
---|
173 | STAT=ierr & |
---|
174 | ) |
---|
175 | IF (ierr /= 0) THEN |
---|
176 | WRITE(*,*) "Could not allocate observation arrays:" |
---|
177 | WRITE(*,*) "dates(", nobs, ")" |
---|
178 | WRITE(*,*) "times(", nobs, ")" |
---|
179 | WRITE(*,*) "lats(", nobs, ")" |
---|
180 | WRITE(*,*) "lons(", nobs, ")" |
---|
181 | WRITE(*,*) "deps(", nlev, ")" |
---|
182 | WRITE(*,*) "variable(", nvar, ")" |
---|
183 | WRITE(*,*) "variable_longname(", nvar, ")" |
---|
184 | WRITE(*,*) "variable_units(", nvar, ")" |
---|
185 | WRITE(*,*) "vals(", nlev, nobs, nvar, ")" |
---|
186 | WRITE(*,*) "model_vals(", nlev, nobs, nvar, ")" |
---|
187 | |
---|
188 | CALL abort |
---|
189 | END IF |
---|
190 | |
---|
191 | |
---|
192 | !!!! set up default values and read in namelist namfbgen.in !!!!!!!!!!!!!!!!!!!! |
---|
193 | ntyp=-1 |
---|
194 | lats(:) = FillValue_real |
---|
195 | lons(:) = FillValue_real |
---|
196 | deps(:) = FillValue_real |
---|
197 | times(:) = FillValue_int |
---|
198 | dates(:) = FillValue_int |
---|
199 | vals(:,:,:) = FillValue_real |
---|
200 | model_vals(:,:,:) = FillValue_real |
---|
201 | variable=REPEAT('X',ilenwmo) |
---|
202 | variable_longname=REPEAT('X',ilenwmo) |
---|
203 | variable_units=REPEAT('X',ilenwmo) |
---|
204 | cdwmo=REPEAT('X',ilenwmo) !station identifier |
---|
205 | cdtyp="90" !station type |
---|
206 | |
---|
207 | WRITE(*,*) "Creating a fdbk file with", nobs, "observations" |
---|
208 | WRITE(*,*) "with", nvar, "variables and",nlev, "depths." |
---|
209 | |
---|
210 | OPEN(10,file='namfbgen.in') |
---|
211 | READ(10,namfbgenerate) |
---|
212 | CLOSE(10) |
---|
213 | |
---|
214 | ! Use bounds to construct the full arrays of coords and values |
---|
215 | |
---|
216 | IF (ln_physpace) THEN |
---|
217 | CALL set_spatial_coords_physpace(lats,lons,nobs,FillValue_real,phys_spacing) |
---|
218 | ELSE IF (ln_gridobs) THEN |
---|
219 | CALL set_spatial_coords_grid(lats,lons,nobs,nlats,nlons,FillValue_real) |
---|
220 | ELSE |
---|
221 | CALL set_spatial_coords(lats,lons,nobs,FillValue_real) |
---|
222 | END IF |
---|
223 | |
---|
224 | ! Add random perturbation to position if ln_perturb_posn=.TRUE. |
---|
225 | IF (ln_perturb_posn) CALL perturb_positions(lats,lons,perturb_limit) |
---|
226 | |
---|
227 | CALL set_datetime(dates,times,julian_dates,nobs,FillValue_int) |
---|
228 | CALL set_depths(deps,nlev,FillValue_real) |
---|
229 | CALL set_obs_values(vals,nvar,nobs,nlev,FillValue_real) |
---|
230 | CALL set_obs_values(model_vals,nvar,nobs,nlev,FillValue_real) |
---|
231 | |
---|
232 | |
---|
233 | !nadd=1 ! Hx |
---|
234 | next=0 ! e.g. TEMP |
---|
235 | ! Add TEMP as extra variable if POTM is defined. |
---|
236 | DO i=1,nvar |
---|
237 | IF (variable(i)(1:4) == "POTM") THEN |
---|
238 | next=1 ! e.g. If set to zero, will not add extra TEMP variable |
---|
239 | EXIT |
---|
240 | END IF |
---|
241 | END DO |
---|
242 | |
---|
243 | CALL init_obfbdata(fbobsdata) |
---|
244 | |
---|
245 | !!!! Allocate the obfb type !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
246 | |
---|
247 | CALL alloc_obfbdata( fbobsdata, nvar, nobs, nlev, nadd, next, .FALSE.) |
---|
248 | |
---|
249 | !!!! Read data specified in the namelist into obfd type!!!!!!!!!!!!!!! |
---|
250 | |
---|
251 | fbobsdata%cname = variable(:) |
---|
252 | fbobsdata%coblong = variable_longname(:) |
---|
253 | fbobsdata%cobunit = variable_units(:) |
---|
254 | fbobsdata%cdjuldref = "19500101000000" |
---|
255 | fbobsdata%cdwmo = cdwmo |
---|
256 | fbobsdata%cdtyp = cdtyp |
---|
257 | |
---|
258 | ! Add model counterpart values if nadd=1 |
---|
259 | IF (nadd > 0) THEN |
---|
260 | IF (nadd /= 1) THEN |
---|
261 | CALL abort |
---|
262 | WRITE(*,*) "Not set-up to add more than 1 additional variable." |
---|
263 | ELSE |
---|
264 | fbobsdata%caddname = "Hx" |
---|
265 | fbobsdata%caddlong(1,:) = "Model interpolated " // variable_longname(:) |
---|
266 | fbobsdata%caddunit(1,:) = variable_units(:) |
---|
267 | fbobsdata%padd(:,:,1,:)=model_vals(:,:,:) |
---|
268 | END IF |
---|
269 | END IF |
---|
270 | |
---|
271 | ! Add TEMP as extra variable if POTM is defined. |
---|
272 | IF (next==1) THEN |
---|
273 | fbobsdata%cextname = "TEMP" |
---|
274 | END IF |
---|
275 | |
---|
276 | ! set up pos/time variables |
---|
277 | fbobsdata%plam(:)=lons(:) |
---|
278 | fbobsdata%pphi(:)=lats(:) |
---|
279 | DO i=1,nobs |
---|
280 | fbobsdata%pdep(:,i)=deps(:) |
---|
281 | END DO |
---|
282 | |
---|
283 | |
---|
284 | fbobsdata%ptim(:)= julian_dates(:) |
---|
285 | |
---|
286 | ! Shuffle the time array to spread the profiles across the time period |
---|
287 | IF (ln_shuffletimes) CALL shuffle(fbobsdata%ptim) |
---|
288 | |
---|
289 | ! read in variable data and flags |
---|
290 | fbobsdata%ivqc(:,:)=1 |
---|
291 | fbobsdata%ivqcf(:,:,:)=1 |
---|
292 | fbobsdata%ivlqc(:,:,:)=1 |
---|
293 | fbobsdata%pob(:,:,:)=vals(:,:,:) |
---|
294 | |
---|
295 | !set up all QC flags |
---|
296 | fbobsdata%ioqc(:)=1 |
---|
297 | fbobsdata%ioqcf(:,:)=0 |
---|
298 | fbobsdata%ipqc(:)=1 |
---|
299 | fbobsdata%ipqcf(:,:)=0 |
---|
300 | fbobsdata%idqc(:,:)=1 !0 |
---|
301 | fbobsdata%idqcf(:,:,:)=0 |
---|
302 | fbobsdata%itqc(:)=1 |
---|
303 | fbobsdata%kindex(:)=0 |
---|
304 | |
---|
305 | !!!! Write out the obfb type !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! |
---|
306 | |
---|
307 | CALL write_obfbdata(TRIM(filenameout),fbobsdata) |
---|
308 | |
---|
309 | DEALLOCATE(dates, & |
---|
310 | times, & |
---|
311 | julian_dates,& |
---|
312 | lats, & |
---|
313 | lons, & |
---|
314 | deps, & |
---|
315 | variable, & |
---|
316 | variable_longname, & |
---|
317 | variable_units, & |
---|
318 | vals, & |
---|
319 | model_vals & |
---|
320 | ) |
---|
321 | |
---|
322 | END IF !test prog |
---|
323 | |
---|
324 | |
---|
325 | END PROGRAM fbgenerate |
---|
326 | |
---|
327 | |
---|