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/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/TOOLS/OBSTOOLS/src – NEMO

source: branches/UKMO/r5518_INGV1_WAVE-coupling/NEMOGCM/TOOLS/OBSTOOLS/src/fbmatchup.F90 @ 7152

Last change on this file since 7152 was 7152, checked in by jcastill, 7 years ago

Initial implementation of wave coupling branch - INGV wave branch + UKMO wave coupling branch

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