New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
fbmatchup.F90 in branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS – NEMO

source: branches/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS/fbmatchup.F90 @ 2893

Last change on this file since 2893 was 2893, checked in by djlea, 13 years ago

Adding obs tools to branch
File size: 14.6 KB
Line 
1PROGRAM fbmatchup
2   USE toolspar_kind
3   USE obs_fbm
4   USE index_sort
5   IMPLICIT NONE
6   !
7   ! Command line arguments for output file and input file
8   !
9#ifndef NOIARGCPROTO
10   INTEGER,EXTERNAL :: iargc
11#endif
12   INTEGER :: nargs
13   CHARACTER(len=256) :: cdoutfile
14   CHARACTER(len=256),ALLOCATABLE :: cdinfile(:)
15   CHARACTER(len=ilenname),ALLOCATABLE :: cdnames(:)
16   INTEGER :: nqc
17   LOGICAL :: ldaily820
18   NAMELIST/namfbmatchup/nqc,ldaily820
19   !
20   ! Input data
21   !
22   TYPE(obfbdata)         :: obsdatatmp(1)
23   TYPE(obfbdata),POINTER :: obsdata(:)
24   INTEGER :: ninfiles,ntotobs,nlev
25   !
26   ! Time sorting arrays
27   !
28   REAL(KIND=dp),ALLOCATABLE :: zsort(:,:)
29   INTEGER,ALLOCATABLE  :: iset(:),inum(:),iindex(:)
30   !
31   ! Comparison arrays and scalars
32   !
33   REAL(KIND=fbsp), ALLOCATABLE :: zrtim(:),zrphi(:),zrlam(:)
34   INTEGER(KIND=SELECTED_INT_KIND(12)), ALLOCATABLE :: irwmo(:)
35   REAL(KIND=fbsp) :: ztim,zphi,zlam
36   INTEGER(KIND=SELECTED_INT_KIND(12)) :: iwmo
37   !
38   ! Output data
39   !
40   TYPE(obfbdata) :: obsoutdata
41   !
42   ! Loop variables
43   !
44   INTEGER :: ia,ip,i1,ii,ij,ik1,ik2,iv,ist
45   LOGICAL :: llfound
46   LOGICAL :: lexists
47   INTEGER :: ityp
48   !
49   ! Get number of command line arguments
50   !
51   nargs = IARGC()
52   IF ( ( MOD(nargs,2) /= 1 ) .OR. ( nargs == 0 )  ) THEN
53      WRITE(*,'(A)')'Usage:'
54      WRITE(*,'(A)')'fbmatchup outputfile inputfile1 varname1 inputfile2 varname2 ...'
55      CALL abort()
56   ENDIF
57   CALL getarg( 1, cdoutfile )
58   !
59   ! Read namelist if present
60   !
61   nqc=1
62   ldaily820=.FALSE.
63   INQUIRE(file='namfbmatchup.in',exist=lexists)
64   IF (lexists) THEN
65      OPEN(10,file='namfbmatchup.in')
66      READ(10,namfbmatchup)
67      CLOSE(10)
68      WRITE(*,namfbmatchup)
69   ENDIF
70   !
71   ! Get input data
72   !
73   ninfiles = ( nargs -1 )/ 2
74   ALLOCATE( obsdata( ninfiles ) )
75   ALLOCATE( cdinfile( ninfiles ) )
76   ALLOCATE( cdnames( ninfiles ) )
77   ip = 1
78   DO ia=1, ninfiles
79      !
80      ! Read the unsorted file
81      !
82      ip = ip + 1
83      CALL getarg( ip, cdinfile(ia) )
84      ip = ip + 1
85      CALL getarg( ip, cdnames(ia) )
86      CALL init_obfbdata( obsdatatmp(1) )
87      CALL read_obfbdata( TRIM(cdinfile(ia)), obsdatatmp(1) )
88      !
89      ! Check that only one variable present in input file
90      !
91      IF ( obsdatatmp(1)%nadd > 1 ) THEN
92         WRITE(*,*)'Warning. More than one variable in input file'
93         WRITE(*,*)'Number of variables = ', obsdatatmp(1)%nadd
94         WRITE(*,*)'Only the first one is going to be stored'
95      ENDIF
96      IF ( obsdatatmp(1)%nadd < 1 ) THEN
97         WRITE(*,*)'Error. less than one variable in input file'
98         CALL abort
99      ENDIF
100      !
101      ! Check that we have few levels than in the first file
102      !
103      IF ( ia > 1 )  THEN
104         IF ( obsdatatmp(1)%nlev > obsdata(1)%nlev ) THEN
105            WRITE(*,*)'Warning. More levels in file than the first file'
106            WRITE(*,*)'Number of levels in current file = ', obsdatatmp(1)%nlev
107            WRITE(*,*)'Number of levels in first file   = ', obsdata(1)%nlev
108            WRITE(*,*)'Only the number of levels in the first'//&
109               &' file will be used'
110         ENDIF
111      ENDIF
112      !
113      ! Check that we have few observations than in the first file
114      !
115      IF ( ia > 1 )  THEN
116         IF ( obsdatatmp(1)%nobs > obsdata(1)%nobs ) THEN
117            WRITE(*,*)'Warning. More obs in file than the first file'
118            WRITE(*,*)'Number of obs in current file = ', obsdatatmp(1)%nobs
119            WRITE(*,*)'Number of obs in first file   = ', obsdata(1)%nobs
120            WRITE(*,*)'Only the observations in the first'//&
121               &' file will be stored'
122         ENDIF
123      ENDIF
124      !
125      ! Check that we have the same number of variables
126      !
127      IF ( ia > 1 )  THEN
128         IF ( obsdatatmp(1)%nvar /= obsdata(1)%nvar ) THEN
129            WRITE(*,*)'Error. Different number of variables.'
130            WRITE(*,*)'Number of var in current file = ', obsdatatmp(1)%nvar
131            WRITE(*,*)'Number of var in first file   = ', obsdata(1)%nvar
132            CALL abort
133         ENDIF
134      ENDIF
135      !
136      ! Check reference datas
137      !
138      IF ( ia > 1 )  THEN
139         IF ( obsdatatmp(1)%cdjuldref /= obsdata(1)%cdjuldref ) THEN
140            WRITE(*,*)'Different reference dates'
141            CALL abort
142         ENDIF
143      ENDIF
144      !
145      ! Special fix for daily average MRB data (820) for the first file
146      !
147      IF (ldaily820.AND.(ia==1)) THEN
148         DO ij = 1,obsdatatmp(1)%nobs
149            READ(obsdatatmp(1)%cdtyp(ij),'(I5)')ityp
150            IF (ityp==820) THEN
151               obsdatatmp(1)%ptim(ij)=INT(obsdatatmp(1)%ptim(ij))+1.0
152            ENDIF
153         ENDDO
154      ENDIF
155      !
156      ! Construct sorting arrays
157      !
158      ALLOCATE( zsort(3,obsdatatmp(1)%nobs), iset(obsdatatmp(1)%nobs), &
159         & inum(obsdatatmp(1)%nobs), iindex(obsdatatmp(1)%nobs))
160      ii = 0
161      DO ij = 1,obsdatatmp(1)%nobs
162         ii = ii+1
163         zsort(1,ii) = obsdatatmp(1)%ptim(ij)
164         zsort(2,ii) = obsdatatmp(1)%pphi(ij)
165         zsort(3,ii) = obsdatatmp(1)%plam(ij)
166         iset(ii) = 1
167         inum(ii) = ij
168      ENDDO
169      !
170      ! Get indexes for time sorting.
171      !
172      CALL index_sort_dp_n(zsort,3,iindex,obsdatatmp(1)%nobs)
173      CALL init_obfbdata( obsdata(ia) )
174      CALL alloc_obfbdata( obsdata(ia), &
175         &                 obsdatatmp(1)%nvar, obsdatatmp(1)%nobs, &
176         &                 obsdatatmp(1)%nlev, obsdatatmp(1)%nadd, &
177         &                 obsdatatmp(1)%next, obsdatatmp(1)%lgrid )
178      !
179      ! Copy input data into output data
180      !
181      CALL merge_obfbdata( 1, obsdatatmp, obsdata(ia), iset, inum, iindex )
182      CALL dealloc_obfbdata( obsdatatmp(1) )
183     
184      WRITE(*,'(2A)')'File = ', TRIM(cdinfile(ia))
185      WRITE(*,'(A,I9,A)')'has', obsdata(ia)%nobs, ' observations'
186
187      DEALLOCATE( zsort, iset, inum, iindex )
188
189   ENDDO
190   !
191   ! Prepare output data
192   !
193   CALL init_obfbdata( obsoutdata )
194   !
195   ! Copy the first input data to output data
196   !
197   CALL copy_obfbdata( obsdata(1), obsoutdata, &
198      &                kadd = ninfiles )
199   DO ia = 1, ninfiles
200      obsoutdata%caddname(ia) = cdnames(ia)
201      obsoutdata%caddlong(ia,:) = obsdata(ia)%caddlong(1,:)
202      obsoutdata%caddunit(ia,:) = obsdata(ia)%caddunit(1,:)
203   ENDDO
204   !
205   ! Allocate comparison arrays and file them
206   !
207   IF (ilenwmo>8) THEN
208      WRITE(*,*)'Fix fbmatchup to allow string length > 8'
209      CALL abort
210   ENDIF
211   ALLOCATE(zrtim(obsoutdata%nobs),zrphi(obsoutdata%nobs), &
212      &     zrlam(obsoutdata%nobs),irwmo(obsoutdata%nobs))
213   DO i1 = 1, obsoutdata%nobs
214      irwmo(i1) = TRANSFER( obsoutdata%cdwmo(i1), irwmo(i1) )
215      zrtim(i1) = REAL( obsoutdata%ptim(i1), fbsp )
216      zrphi(i1) = REAL( obsoutdata%pphi(i1), fbsp )
217      zrlam(i1) = REAL( obsoutdata%plam(i1), fbsp ) 
218   ENDDO
219   !
220   ! Merge extra data into output data
221   !
222   DO ia = 2, ninfiles
223      ist = 1
224      DO ii = 1, obsdata(ia)%nobs
225         IF (MOD(ii,10000)==1) THEN
226            WRITE(*,*)'Handling observation no ',ii,' for file no ',ia
227         ENDIF
228         llfound = .FALSE.
229         iwmo = TRANSFER( obsdata(ia)%cdwmo(ii), iwmo )
230         ztim = REAL( obsdata(ia)%ptim(ii), fbsp )
231         zphi = REAL( obsdata(ia)%pphi(ii), fbsp )
232         zlam = REAL( obsdata(ia)%plam(ii), fbsp ) 
233         ! Check if the the same index is the right one.
234         IF ( iwmo == irwmo(ii) ) THEN
235            IF ( ztim == zrtim(ii) ) THEN
236               IF (  zphi == zrphi(ii) ) THEN
237                  IF ( zlam == zrlam(ii) ) THEN
238                     llfound = .TRUE.
239                     DO iv = 1, obsdata(ia)%nvar
240                        ! Since the inner loop don't change the
241                        ! qc decisions use this to ensure match
242                        ! for duplicate observations
243                        IF ( obsdata(ia)%ivqc(ii,iv) /= &
244                           & obsoutdata%ivqc(ii,iv)) THEN
245                           llfound = .FALSE.
246                           CYCLE
247                        ENDIF
248                     ENDDO
249                     IF (llfound) i1 = ii
250                  ENDIF
251               ENDIF
252            ENDIF
253         ENDIF
254         ! Search for position in from previous found position
255         ! if not the same index
256         IF (.NOT.llfound) THEN
257            DO i1 = ist, obsoutdata%nobs
258               IF ( iwmo == irwmo(i1) ) THEN
259                  IF ( ztim == zrtim(i1) ) THEN
260                     IF (  zphi == zrphi(i1) ) THEN
261                        IF ( zlam == zrlam(i1) ) THEN
262                           llfound = .TRUE.
263                           DO iv = 1, obsdata(ia)%nvar
264                              ! Since the inner loop don't change the
265                              ! qc decisions use this to ensure match
266                              ! for duplicate observations
267                              IF ( obsdata(ia)%ivqc(ii,iv) /= &
268                                 & obsoutdata%ivqc(i1,iv)) THEN
269                                 llfound = .FALSE.
270                                 WRITE (*,*)'QC flags different for ',&
271                                    &  TRIM(obsdata(ia)%cdwmo(ii)),' at ', &
272                                    & obsdata(ia)%ptim(ii)
273                                 CYCLE
274                              ENDIF
275                           ENDDO
276                           IF (llfound) EXIT
277                        ENDIF
278                     ENDIF
279                  ENDIF
280               ENDIF
281            ENDDO
282         ENDIF
283         ! If not fount try agan from the beginnning
284         IF ( .NOT.llfound ) THEN
285            DO i1 = 1, obsoutdata%nobs
286               IF ( iwmo == irwmo(i1) ) THEN
287                  IF ( ztim == zrtim(i1) ) THEN
288                     IF (  zphi == zrphi(i1) ) THEN
289                        IF ( zlam == zrlam(i1) ) THEN
290                           llfound = .TRUE.
291                           DO iv = 1, obsdata(ia)%nvar
292                              ! Since the inner loop don't change the
293                              ! qc decisions use this to ensure match
294                              ! for duplicate observations
295                              IF ( obsdata(ia)%ivqc(ii,iv) /= &
296                                 & obsoutdata%ivqc(i1,iv)) THEN
297                                 llfound = .FALSE.
298                                 WRITE (*,*)'QC flags different for ',&
299                                    &  TRIM(obsdata(ia)%cdwmo(ii)),' at ', &
300                                    & obsdata(ia)%ptim(ii)
301                                 CYCLE
302                              ENDIF
303                           ENDDO
304                           IF (llfound) EXIT
305                        ENDIF
306                     ENDIF
307                  ENDIF
308               ENDIF
309            ENDDO
310         ENDIF
311         ! If found put the data into the common structure
312         IF (llfound) THEN
313            IF ( nqc == ia ) THEN
314               obsoutdata%ioqc(i1)   = obsdata(ia)%ioqc(ii)
315               obsoutdata%ipqc(i1)   = obsdata(ia)%ipqc(ii)
316               obsoutdata%itqc(i1)   = obsdata(ia)%itqc(ii)
317               obsoutdata%ivqc(i1,:) = obsdata(ia)%ivqc(ii,:)
318            ENDIF
319            obsoutdata%ioqcf(:,i1)   = IOR( obsdata(ia)%ioqcf(:,ii), &
320               &                            obsoutdata%ioqcf(:,i1) )
321            obsoutdata%ipqcf(:,i1)   = IOR( obsdata(ia)%ipqcf(:,ii), &
322               &                            obsoutdata%ipqcf(:,i1) )
323            obsoutdata%itqcf(:,i1)   = IOR( obsdata(ia)%itqcf(:,ii), &
324               &                            obsoutdata%itqcf(:,i1) )
325            obsoutdata%ivqcf(:,i1,:) = IOR( obsdata(ia)%ivqcf(:,ii,:), &
326               &                            obsoutdata%ivqcf(:,i1,:) )
327            llfound = .FALSE.
328            ! Search for levels
329            DO ik1 = 1, obsdata(ia)%nlev
330               DO ik2 = 1, obsoutdata%nlev
331                  IF ( REAL( obsdata(ia)%pdep(ik1,ii), fbsp ) ==  &
332                     & REAL( obsoutdata%pdep(ik2,i1),  fbsp ) ) THEN
333                     obsoutdata%padd(ik2,i1,ia,:) = &
334                        & obsdata(ia)%padd(ik1,ii,1,:)
335                     IF ( nqc == ia ) THEN
336                        obsoutdata%idqc(ik2,i1)    = obsdata(ia)%idqc(ik1,ii)
337                        obsoutdata%ivlqc(ik2,i1,:) = obsdata(ia)%ivlqc(ik1,ii,:)
338                     ENDIF
339                     obsoutdata%idqcf(:,ik2,i1)    = &
340                        &                IOR( obsdata(ia)%idqcf(:,ik1,ii), &
341                        &                     obsoutdata%idqcf(:,ik2,i1) )
342                     obsoutdata%ivlqcf(:,ik2,i1,:) = &
343                        &                IOR( obsdata(ia)%ivlqcf(:,ik1,ii,:), &
344                        &                     obsoutdata%ivlqcf(:,ik2,i1,:) )
345                     llfound = .TRUE.
346                     EXIT
347                  ENDIF
348               ENDDO
349               ! Write warning if level not found
350               IF (.NOT.llfound.AND.(obsdata(ia)%pdep(ik1,ii)/=fbrmdi)) THEN
351                  WRITE(*,*)'Level not found in first file : ',&
352                     &      TRIM( cdinfile(1)  )
353                  WRITE(*,*)'Data file                     : ',&
354                     &      TRIM( cdinfile(ia) )
355                  WRITE(*,*)'Identifier                    : ',&
356                     &      obsdata(ia)%cdwmo(ii)
357                  WRITE(*,*)'Julian date                   : ',&
358                     &      obsdata(ia)%ptim(ii)
359                  WRITE(*,*)'Latitude                      : ',&
360                     &      obsdata(ia)%pphi(ii)
361                  WRITE(*,*)'Longitude                     : ',&
362                     &      obsdata(ia)%plam(ii)
363                  WRITE(*,*)'Depth                         : ',&
364                     &      obsdata(ia)%pdep(ik1,ii)
365               ENDIF
366            ENDDO
367            ist = i1
368         ELSE
369            ! Write warning if observation not found
370            WRITE(*,*)'Observation not found in first data file : ',&
371               &      TRIM( cdinfile(1)  )
372            WRITE(*,*)'Data file                                : ',&
373               &      TRIM( cdinfile(ia) )
374            WRITE(*,*)'Identifier                               : ',&
375               &      obsdata(ia)%cdwmo(ii)
376            WRITE(*,*)'Julian date                              : ',&
377               &      obsdata(ia)%ptim(ii)
378            WRITE(*,*)'Latitude                                 : ',&
379               &      obsdata(ia)%pphi(ii)
380            WRITE(*,*)'Longitude                                : ',&
381               &      obsdata(ia)%plam(ii)
382            ist = 1
383         ENDIF
384      ENDDO
385      IF (obsdata(ia)%nobs>0) THEN
386         WRITE(*,*)'Handled last obs. no    ',ii,' for file no ',ia
387      ENDIF
388   ENDDO
389   !
390   ! Write output file
391   !
392   CALL write_obfbdata( TRIM(cdoutfile), obsoutdata )
393
394END PROGRAM fbmatchup
Note: See TracBrowser for help on using the repository browser.