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.
fbsel.F90 in branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/TOOLS/OBSTOOLS/src – NEMO

source: branches/2014/dev_r4650_UKMO7_STARTHOUR/NEMOGCM/TOOLS/OBSTOOLS/src/fbsel.F90 @ 5985

Last change on this file since 5985 was 5985, checked in by timgraham, 8 years ago

Reinstate keywords before upgrading to head of trunk
  • Property svn:keywords set to Id
File size: 8.7 KB
Line 
1PROGRAM fbsel
2   !!---------------------------------------------------------------------
3   !!
4   !!                     ** PROGRAM fbsel **
5   !!
6   !!  ** Purpose : Select or subsample observations
7   !!
8   !!  ** Method  : Use of utilities from obs_fbm.
9   !!
10   !!  ** Action  :
11   !!
12   !!   Usage:
13   !!     fbsel.exe <input filename> <output filename>
14   !!
15   !!   History :
16   !!        ! 2010 (K. Mogensen) Initial version
17   !!----------------------------------------------------------------------
18   USE obs_fbm
19   USE date_utils
20   IMPLICIT NONE
21   TYPE(obfbdata) :: fbdatain,fbdataout
22   CHARACTER(len=256) :: filenamein,filenameout,filenametmp,cnameout
23#ifndef NOIARGCPROTO
24   INTEGER,EXTERNAL :: iargc
25#endif
26   INTEGER,PARAMETER :: maxtyp=1023
27   INTEGER,PARAMETER :: maxdates=20
28   INTEGER :: nqc,ntyp,ndates,ninidate(maxdates),nenddate(maxdates)
29   LOGICAL :: lsplitqc,lsplittyp,lsplitstat
30   INTEGER :: iqc,ityp,idate,istat
31   REAL :: maxlat,minlat,maxlon,minlon
32   CHARACTER(len=ilenwmo) :: cdwmo,cdwmobeg,cdwmoend
33   CHARACTER(len=ilenwmo), DIMENSION(:), POINTER :: clstatids
34   INTEGER :: nstat
35   NAMELIST/namsel/nqc,ntyp,ndates,ninidate,nenddate,lsplitqc,lsplittyp, &
36      &            lsplitstat,maxlat,minlat,maxlon,minlon,cdwmo,&
37      &            cdwmobeg,cdwmoend
38
39   IF (iargc()/=2) THEN
40      WRITE(*,*)'Usage:'
41      WRITE(*,*)'fbsel <input filename> <output filename>'
42      CALL abort
43   ENDIF
44
45   CALL getarg(1,filenamein)
46   CALL getarg(2,filenameout)
47   
48   nqc=-1
49   ntyp=-1
50   ndates=1
51   ninidate=19500101
52   nenddate=21000101
53
54   lsplitqc=.FALSE.
55   lsplittyp=.FALSE.
56   lsplitstat=.FALSE.
57   cdwmo=REPEAT('X',ilenwmo)
58   cdwmobeg=cdwmo
59   cdwmoend=cdwmo
60   maxlat=1e+38
61   minlat=-1e+38
62   maxlon=1e+38
63   minlon=-1e+38
64   OPEN(10,file='namsel.in')
65   READ(10,namsel)
66   CLOSE(10)
67   IF (cdwmobeg==REPEAT('X',ilenwmo)) cdwmobeg=cdwmo
68   IF (cdwmoend==REPEAT('X',ilenwmo)) cdwmoend=cdwmo
69   WRITE(*,namsel)
70
71   CALL init_obfbdata(fbdatain)
72   CALL init_obfbdata(fbdataout)
73
74   WRITE(*,*)'Reading file : ',TRIM(filenamein)
75   CALL read_obfbdata(TRIM(filenamein),fbdatain)
76   WRITE(*,*)'Number of observations before selection = ',fbdatain%nobs
77   DO idate=1,ndates
78      IF (ndates==1) THEN
79         cnameout=filenameout
80      ELSE
81         WRITE(cnameout,'(I2.2,2A)')idate,'_',TRIM(filenameout)
82      ENDIF
83      IF (lsplitqc) THEN
84         DO iqc=1,3
85            CALL fb_sel(fbdatain,fbdataout,iqc,ntyp, &
86               &        ninidate(idate),nenddate(idate), &
87               &        maxlat,minlat,maxlon,minlon,cdwmobeg,cdwmoend)
88            WRITE(filenametmp,'(A,I2.2,A,A)')'qc_',iqc,'_',TRIM(cnameout)
89            IF (fbdataout%nobs>0) THEN
90               WRITE(*,*)'QC selected = ',iqc
91               WRITE(*,*)'Number of observations selected = ',fbdataout%nobs
92               WRITE(*,*)'Writing file : ',TRIM(filenametmp)
93               CALL write_obfbdata(TRIM(filenametmp),fbdataout)
94            ENDIF
95            CALL dealloc_obfbdata(fbdataout)
96         ENDDO
97      ELSEIF (lsplittyp) THEN
98         DO ityp=0,maxtyp
99            CALL fb_sel(fbdatain,fbdataout,nqc,ityp, &
100               &        ninidate(idate),nenddate(idate), &
101               &        maxlat,minlat,maxlon,minlon,cdwmobeg,cdwmoend)
102            WRITE(filenametmp,'(A,I4.4,A,A)')'typ_',ityp,'_',TRIM(cnameout)
103            IF (fbdataout%nobs>0) THEN
104               WRITE(*,*)'Type = ',ityp
105               WRITE(*,*)'Number of observations selected = ',fbdataout%nobs
106               WRITE(*,*)'Writing file : ',TRIM(filenametmp)
107               CALL write_obfbdata(TRIM(filenametmp),fbdataout)
108            ENDIF
109            CALL dealloc_obfbdata(fbdataout)
110         ENDDO
111      ELSEIF (lsplitstat) THEN
112         CALL fb_sel_uniqueids(fbdatain,clstatids,nstat)
113         DO istat=1,nstat
114            CALL fb_sel(fbdatain,fbdataout,nqc,ntyp, &
115               &        ninidate(idate),nenddate(idate), &
116               &        maxlat,minlat,maxlon,minlon,clstatids(istat),clstatids(istat))
117            WRITE(filenametmp,'(4A)')'statid_', &
118               & TRIM(clstatids(istat)),'_',TRIM(cnameout)
119            IF (fbdataout%nobs>0) THEN
120               WRITE(*,*)'Station = ',clstatids(istat)
121               WRITE(*,*)'Number of observations selected = ',fbdataout%nobs
122               WRITE(*,*)'Writing file : ',TRIM(filenametmp)
123               CALL write_obfbdata(TRIM(filenametmp),fbdataout)
124            ENDIF
125            CALL dealloc_obfbdata(fbdataout)
126         ENDDO
127      ELSE
128         CALL fb_sel(fbdatain,fbdataout,nqc,ntyp, & 
129            &        ninidate(idate),nenddate(idate), &
130            &        maxlat,minlat,maxlon,minlon,cdwmobeg,cdwmoend)
131         WRITE(*,*)'Number of observations selected = ',fbdataout%nobs
132         WRITE(*,*)'Writing file : ',TRIM(cnameout)
133         CALL write_obfbdata(TRIM(cnameout),fbdataout)
134         CALL dealloc_obfbdata(fbdataout)
135      ENDIF
136   ENDDO
137
138CONTAINS
139
140   SUBROUTINE fb_sel(fbdatain,fbdataout,nqc,ntyp,ninidate,nenddate,&
141      &              maxlat,minlat,maxlon,minlon,cdwmobeg,cdwmoend)
142      TYPE(obfbdata) :: fbdatain,fbdataout
143      INTEGER :: nqc,ntyp,ninidate,nenddate
144      REAL :: maxlat,minlat,maxlon,minlon
145      CHARACTER(len=ilenwmo) :: cdwmobeg,cdwmoend
146      INTEGER :: jobs
147      INTEGER :: iqc,ityp
148      LOGICAL :: llvalid(fbdatain%nobs)
149      INTEGER :: iyea,imon,iday
150      REAL(KIND=dp) :: zjini,zjend
151      LOGICAL :: lcheckwmo
152
153      lcheckwmo=(cdwmobeg/=REPEAT('X',ilenwmo)).OR.&
154         &      (cdwmoend/=REPEAT('X',ilenwmo))
155      iyea=ninidate/10000
156      imon=ninidate/100-iyea*100
157      iday=ninidate-iyea*10000-imon*100
158      CALL greg2jul(0,0,0,iday,imon,iyea,zjini)
159      iyea=nenddate/10000
160      imon=nenddate/100-iyea*100
161      iday=nenddate-iyea*10000-imon*100
162      CALL greg2jul(0,0,0,iday,imon,iyea,zjend)
163      DO jobs = 1, fbdatain%nobs
164         llvalid(jobs)=.TRUE.
165         IF (nqc/=-1) THEN
166            CALL check_prof(fbdatain,jobs,iqc)
167            llvalid(jobs)=(iqc==nqc).AND.llvalid(jobs)
168         ENDIF
169         IF (ntyp/=-1) THEN
170            READ(fbdatain%cdtyp(jobs),'(I4)')ityp
171            llvalid(jobs)=(ityp==ntyp).AND.llvalid(jobs)
172         ENDIF
173         IF (ninidate/=-1) THEN
174            llvalid(jobs)=(fbdatain%ptim(jobs)>zjini).AND.llvalid(jobs)
175         ENDIF
176         IF (nenddate/=-1) THEN
177            llvalid(jobs)=(fbdatain%ptim(jobs)<=zjend).AND.llvalid(jobs)
178         ENDIF
179         llvalid(jobs)=(fbdatain%pphi(jobs)<=maxlat).AND.       &
180            &          (fbdatain%pphi(jobs)>=minlat).AND.       &
181            &          (((fbdatain%plam(jobs)<=maxlon).AND.     &
182            &            (fbdatain%plam(jobs)>=minlon)).OR.     &
183            &           ((fbdatain%plam(jobs)+360<=maxlon).AND. &
184            &            (fbdatain%plam(jobs)+360>=minlon)).OR. &
185            &           ((fbdatain%plam(jobs)-360<=maxlon).AND. &
186            &            (fbdatain%plam(jobs)-360>=minlon))).AND.llvalid(jobs) 
187         IF (lcheckwmo) THEN
188            llvalid(jobs)=LGE(TRIM(fbdatain%cdwmo(jobs)),TRIM(cdwmobeg)) &
189               &    .AND. LLE(TRIM(fbdatain%cdwmo(jobs)),TRIM(cdwmoend)) &
190               &    .AND. llvalid(jobs)
191         ENDIF
192         ! Add more checks here...
193      ENDDO
194
195      CALL subsamp_obfbdata(fbdatain,fbdataout,llvalid)
196
197   END SUBROUTINE fb_sel
198
199   SUBROUTINE fb_sel_uniqueids(fbdatain,clstatids,nstat)
200      TYPE(obfbdata) :: fbdatain
201      CHARACTER(len=ilenwmo), DIMENSION(:), POINTER :: clstatids
202      INTEGER :: nstat
203      INTEGER :: jobs,kobs
204      LOGICAL, DIMENSION(fbdatain%nobs) :: lunique
205
206      lunique(:)=.TRUE.
207      DO jobs=1,fbdatain%nobs
208         IF (.NOT.lunique(jobs)) CYCLE
209         DO kobs=jobs+1,fbdatain%nobs
210            IF (.NOT.lunique(kobs)) CYCLE
211            IF (fbdatain%cdwmo(jobs)==fbdatain%cdwmo(kobs)) THEN
212               lunique(kobs)=.FALSE.
213            ENDIF
214         ENDDO
215      ENDDO
216      nstat=COUNT(lunique)
217      ALLOCATE(clstatids(nstat))
218      kobs=0
219      DO jobs=1,fbdatain%nobs
220         IF (lunique(jobs)) THEN
221            kobs=kobs+1
222            clstatids(kobs)=fbdatain%cdwmo(jobs)
223         ENDIF
224      ENDDO
225      WRITE(*,*)'Unique station ids'
226      DO jobs=1,nstat
227         WRITE(*,'(I5,1X,A)')jobs,clstatids(jobs)
228      ENDDO
229
230   END SUBROUTINE fb_sel_uniqueids
231   
232   SUBROUTINE check_prof(fbdata,iobs,iqc)
233     
234      TYPE(obfbdata) :: fbdata
235      INTEGER :: iobs,iqc
236      INTEGER :: i,ivar
237     
238      LOGICAL :: lpart,lfull
239
240      lpart=.false.
241      lfull=.true.
242      DO ivar=1,fbdata%nvar
243         DO i=1,fbdata%nlev
244            IF ((fbdata%ivlqc(i,iobs,ivar)>2).AND.&
245               &(fbdata%ivlqc(i,iobs,ivar)<9)) lpart = .TRUE.
246            IF (fbdata%ivlqc(i,iobs,ivar)<=2) lfull = .FALSE.
247         ENDDO
248      ENDDO
249
250      IF(lfull) THEN
251         iqc=3
252      ELSEIF (lpart) then
253         iqc=2
254      ELSE
255         iqc=1
256      ENDIF
257
258   END SUBROUTINE check_prof
259
260END PROGRAM fbsel
Note: See TracBrowser for help on using the repository browser.