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.
fbstat.F90 in utils/tools/OBSTOOLS/src – NEMO

source: utils/tools/OBSTOOLS/src/fbstat.F90 @ 10841

Last change on this file since 10841 was 3002, checked in by djlea, 12 years ago

Update documentation for obstools and dataplot. Removal of dataplot code not needed. Addition of headers to some dataplot code. Addition of .exe to command example in obstools.

  • Property svn:executable set to *
File size: 51.7 KB
Line 
1PROGRAM fbstat
2   !!---------------------------------------------------------------------
3   !!
4   !!                     ** PROGRAM fbstat **
5   !!
6   !!  ** Purpose : Output feedback file summary info/statistics
7   !!         into a number of .dat files for different areas
8   !!
9   !!  ** Method  : Use of utilities from obs_fbm.
10   !!
11   !!  ** Action  :
12   !!
13   !!   Usage:
14   !!     fbstat.exe [-nmlev] <filenames>
15   !!   Optional:
16   !!     namelist = namfbstat.in
17   !!
18   !!   History :
19   !!        ! 2010 (K. Mogensen) Initial version
20   !!----------------------------------------------------------------------
21   USE obs_fbm
22   USE fbaccdata
23   USE coords
24   USE omonainfo
25   USE fbstatncio
26   USE proftools
27   IMPLICIT NONE
28   TYPE(obfbdata) :: fbdata
29   CHARACTER(len=256) :: filename,outfilename
30   INTEGER :: jfile,jbox,jlev,jfirst,jvar,jadd,ji,ja,jt,jboxl
31#ifndef NOIARGCPROTO
32   INTEGER,EXTERNAL :: iargc
33#endif
34   REAL,DIMENSION(:),ALLOCATABLE :: zlev
35   INTEGER :: nmlev, nfiles
36   LOGICAL :: lexists,lomona,ltext,lnetcdf,lzinp
37   REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zdat3d
38   REAL, ALLOCATABLE, DIMENSION(:,:) :: zdat2d
39   INTEGER,DIMENSION(1) :: itime
40   INTEGER :: inidate,icurdate,loopno,ityp,iloopno
41   INTEGER :: nvar,nadd,noberr,nbgerr
42   CHARACTER(len=4) :: expver
43   CHARACTER(len=20) :: cltyp
44   CHARACTER(len=128) :: cdfmthead,cdfmtbody
45   LOGICAL :: lnear,linner,linnerp,linnerini,lpassive,lhistogram,lfound
46   LOGICAL :: lxyplot,lrmmean
47   INTEGER :: nqc,nqco
48   REAL :: rlspc,rlmax
49   CHARACTER(len=ilenname), DIMENSION(:), ALLOCATABLE :: cname,caddname,&
50      & cobename,cbgename
51   INTEGER, PARAMETER :: nmaxareas = 20
52   CHARACTER(len=20), DIMENSION(nmaxareas) :: carea
53   LOGICAL, DIMENSION(:), ALLOCATABLE :: lskipbox
54   INTEGER, parameter :: maxtyp = 10
55   CHARACTER(len=ilentyp), DIMENSION(maxtyp) :: ctyp
56   INTEGER :: ntyp,nboxuserl,ipdcst
57   REAL :: mindcst
58   NAMELIST/namfbstat/ltext,lomona,lnetcdf,nmlev,inidate,icurdate,loopno,&
59      &               expver,lnear,linner,lpassive,lhistogram,&
60      &               zhistmax,zhistmin,zhiststep,zcheck,carea,nmlev,&
61      &               nqc,nqco, &
62      &               rlspc,rlmax,ntyp,ctyp,&
63      &               lxyplot,zxymin,zxymax,zxystep,lzinp,lrmmean,mindcst
64
65   ltext=.TRUE.
66   lnetcdf=.TRUE.
67   lomona=.FALSE.
68   nmlev=31
69   inidate=19010101
70   icurdate=19010116
71   loopno=0
72   expver='test'
73   lnear=.TRUE.
74   linner=.FALSE.
75   lpassive=.FALSE.
76   lhistogram=.FALSE.
77   zhistmin(:)=-10.0
78   zhistmax(:)=10.0
79   zhiststep(:)=0.1
80   zcheck(:)=1000.0
81   nqc=2
82   nqco=2
83   carea(:)='all'
84   rlmax=5000.0
85   rlspc=-0.1
86   ntyp=1
87   ctyp(:)='all'
88   lxyplot=.FALSE.
89   zxymin(:)=-5.0
90   zxymax(:)=45.0
91   zxystep(:)=0.5
92   lzinp=.FALSE.
93   lrmmean=.FALSE.
94   mindcst=-1.0
95
96   INQUIRE(file='namfbstat.in',exist=lexists)
97   IF (lexists) THEN
98      OPEN(10,file='namfbstat.in')
99      READ(10,namfbstat)
100      CLOSE(10)
101      WRITE(*,namfbstat)
102   ENDIF
103   mindcst=mindcst*1000.0 !From km to m.
104   IF (iargc()==0) THEN
105      WRITE(*,*)'Usage:'
106      WRITE(*,*)'fbstat [-nmlev] <filenames>'
107      CALL abort
108   ENDIF
109   jfirst=1
110   DO ji=1,2
111      CALL getarg(jfirst,filename)
112      IF (filename=='-42') THEN
113         nmlev=42
114         jfirst=jfirst+1
115      ELSEIF(filename=='-31') THEN
116         nmlev=31
117         jfirst=jfirst+1
118      ELSEIF(filename=='-1') THEN
119         nmlev=1
120         lnear=.TRUE.
121         jfirst=jfirst+1
122      ELSEIF(filename=='-q') THEN
123         jfirst=jfirst+1
124         CALL getarg(jfirst,filename)
125         READ(filename,'(I4)')nqc
126         IF ((nqc<0).OR.(nqc>4)) THEN
127            WRITE(*,*)'Quality control option (-q) should be 1 to 4'
128            CALL abort
129         ENDIF
130         jfirst=jfirst+1
131      ENDIF
132   END DO
133   nfiles=iargc()
134
135   CALL coord_user_init('o')
136
137   ALLOCATE(lskipbox(nboxuser))
138   lskipbox(:)=.FALSE.
139
140   IF (carea(1)/='all') THEN
141      IF (lomona) THEN
142         WRITE(*,*)'For omona files carea(1) has to be all'
143         CALL abort
144      ENDIF
145      lskipbox(:)=.TRUE.
146      DO ji=1,nmaxareas
147         IF (carea(ji)/='all') THEN
148            lfound=.FALSE.
149            DO jbox=1,nboxuser
150               IF (TRIM(carea(ji))==TRIM(cl_boxes_user(jbox))) THEN
151                  lskipbox(jbox)=.FALSE.
152                  lfound=.TRUE.
153               ENDIF
154            ENDDO
155            IF (.NOT.lfound) THEN
156               WRITE(*,*)'Area ',TRIM(carea(ji)),' not found'
157               CALL abort
158            ENDIF
159         ENDIF
160      ENDDO
161      nboxuserl=0
162      DO ji=1,nboxuser
163         WRITE(*,*)'Area ',TRIM(cl_boxes_user(ji)),' is set to ',lskipbox(ji)
164         IF (.NOT.lskipbox(ji)) nboxuserl=nboxuserl+1
165      ENDDO
166      WRITE(*,*)'Total areas for statistics = ',nboxuserl
167      IF (lomona.AND.(nboxuserl/=nboxuser)) THEN
168         WRITE(*,*)'Omona files only possible if all areas'
169         CALL abort
170      ENDIF
171   ELSE
172      nboxuserl=nboxuser
173   ENDIF
174
175   IF (rlspc>0.0) THEN
176      lnear=.TRUE.
177      nmlev=rlmax/rlspc+1
178      ALLOCATE(zlev(nmlev))
179      DO ji=1,nmlev
180         zlev(ji)=(ji-1)*rlspc
181      ENDDO
182   ELSE
183      IF (.NOT.lnear) nmlev=nmlev-1
184      ALLOCATE(&
185         & zlev(nmlev) &
186         & )
187      IF(lnear) THEN
188         CALL getlevs(nmlev,zlev)
189      ELSE
190         CALL getlevsmean(nmlev,zlev)
191      ENDIF
192   ENDIF
193
194   DO jfile=jfirst, nfiles
195      CALL getarg(jfile,filename)
196      WRITE(*,*)'Handling file : ',TRIM(filename)
197      CALL flush(6)
198      IF (lzinp) THEN
199#if defined NOSYSTEM
200         WRITE(*,*)'Compressed files need the system subroutine call'
201         CALL abort
202#else
203         CALL system('cp '//TRIM(filename)//' fbstat_tmp.nc.gz')
204         CALL system('gzip -df fbstat_tmp.nc.gz')
205         CALL read_obfbdata('fbstat_tmp.nc',fbdata)
206         CALL system('rm -f fbstat_tmp.nc')
207#endif
208      ELSE
209         CALL read_obfbdata(TRIM(filename),fbdata)
210      ENDIF
211      CALL sealsfromargo( fbdata )
212      IF (jfile==jfirst) THEN
213         nvar=fbdata%nvar
214         nadd=0
215         DO ja= 1, fbdata%nadd
216            IF (fbdata%caddname(ja)(1:2)=='Hx') nadd=nadd+1
217         ENDDO
218         noberr=0
219         DO ja= 1, fbdata%nadd
220            IF (fbdata%caddname(ja)(1:3)=='OBE') noberr=noberr+1
221         ENDDO
222         nbgerr=0
223         DO ja= 1, fbdata%nadd
224            IF (fbdata%caddname(ja)(1:3)=='BGE') nbgerr=nbgerr+1
225         ENDDO
226         IF (lhistogram) THEN
227            IF (nvar>maxvars) THEN
228               WRITE(*,*)'fbstat.F90: Increase maxvars to ',nvar
229               WRITE(*,*)'if you want histograms'
230               CALL abort
231            ENDIF
232            DO jvar = 1, nvar
233               hist(jvar)%npoints=(zhistmax(jvar)-zhistmin(jvar))&
234                  &              /zhiststep(jvar)+1
235               WRITE(*,*)'Number of points in histogram      = ',&
236                  &  hist(jvar)%npoints
237               WRITE(*,*)'Size of histogram array (elements) = ',&
238                  &  hist(jvar)%npoints*nmlev*nboxuserl*nadd
239               ALLOCATE(&
240                  & hist(jvar)%nhist(hist(jvar)%npoints,nmlev,nboxuserl,nadd,ntyp) &
241                  & )
242               hist(jvar)%nhist(:,:,:,:,:)=0
243            ENDDO
244         ENDIF
245         ipdcst=0
246         IF (mindcst>0) THEN
247            DO ja= 1, fbdata%next
248               IF (TRIM(fbdata%cextname(ja))=='DCST') THEN
249                  ipdcst=ja
250                  EXIT
251               ENDIF
252            ENDDO
253            IF (ipdcst==0) THEN
254               WRITE(*,*)'Distance to coast not found in file, but mindcst>0'
255               WRITE(*,*)'Extra variables:'
256               DO ja= 1, fbdata%next
257                  WRITE(*,*)ja,TRIM(fbdata%cextname(ja))
258               ENDDO
259               CALL abort
260            ENDIF
261         ENDIF
262         IF (lxyplot) THEN
263            IF (nvar>maxvars) THEN
264               WRITE(*,*)'fbstat.F90: Increase maxvars to ',nvar
265               WRITE(*,*)'if you want xyplots'
266               CALL abort
267            ENDIF
268            DO jvar = 1, nvar
269               xy(jvar)%npoints=(zxymax(jvar)-zxymin(jvar))&
270                  &              /zxystep(jvar)+1
271               WRITE(*,*)'Number of points in x and y for xyplots = ',&
272                  &  xy(jvar)%npoints
273               WRITE(*,*)'Size of xyplot array (elements)         = ',&
274                  &  xy(jvar)%npoints*xy(jvar)%npoints*nmlev*nboxuserl*nadd
275               ALLOCATE(&
276                  & xy(jvar)%nxy(xy(jvar)%npoints,xy(jvar)%npoints,&
277                  &              nmlev,nboxuserl,nadd,ntyp) &
278                  & )
279               xy(jvar)%nxy(:,:,:,:,:,:)=0
280            ENDDO
281         ENDIF
282         ALLOCATE(&
283            & inum(nmlev,nboxuserl,nadd,nvar,ntyp),  &
284            & inumov(nmlev,nboxuserl,noberr,nvar,ntyp),  &
285            & inumbv(nmlev,nboxuserl,nbgerr,nvar,ntyp),  &
286            & inuma(nmlev,nboxuserl,nvar,ntyp),      &
287            & zbias(nmlev,nboxuserl,nadd,nvar,ntyp), &
288            & zrms(nmlev,nboxuserl,nadd,nvar,ntyp),  &
289            & zsdev(nmlev,nboxuserl,nadd,nvar,ntyp), &
290            & zomean(nmlev,nboxuserl,nadd,nvar,ntyp),&
291            & zmmean(nmlev,nboxuserl,nadd,nvar,ntyp),&
292            & zoemea(nmlev,nboxuserl,noberr,nvar,ntyp),&
293            & zovmea(nmlev,nboxuserl,noberr,nvar,ntyp),&
294            & zbemea(nmlev,nboxuserl,nbgerr,nvar,ntyp),&
295            & zbvmea(nmlev,nboxuserl,nbgerr,nvar,ntyp),&
296            & zoamean(nmlev,nboxuserl,nvar,ntyp),    &
297            & cname(nvar),                           &
298            & caddname(nadd),                        &
299            & cobename(noberr),                      &
300            & cbgename(nbgerr)                       &
301            & )
302         DO jvar = 1, nvar
303            cname(jvar) = fbdata%cname(jvar)
304         END DO
305         jadd = 0
306         DO ja= 1, fbdata%nadd
307            IF (fbdata%caddname(ja)(1:2)=='Hx') THEN
308               jadd=jadd+1
309               caddname(jadd) = fbdata%caddname(ja)
310            ENDIF
311         END DO
312         jadd = 0
313         DO ja= 1, fbdata%nadd
314            IF (fbdata%caddname(ja)(1:3)=='OBE') THEN
315               jadd=jadd+1
316               cobename(jadd) = fbdata%caddname(ja)
317            ENDIF
318         END DO
319         jadd = 0
320         DO ja= 1, fbdata%nadd
321            IF (fbdata%caddname(ja)(1:3)=='BGE') THEN
322               jadd=jadd+1
323               cbgename(jadd) = fbdata%caddname(ja)
324            ENDIF
325         END DO
326         IF (nadd>0) THEN
327            inum(:,:,:,:,:)=0
328            zbias(:,:,:,:,:)=0.0
329            zrms(:,:,:,:,:)=0.0
330            zsdev(:,:,:,:,:)=0.0
331            zomean(:,:,:,:,:)=0.0 
332            zmmean(:,:,:,:,:)=0.0
333         ENDIF
334         IF (noberr>0) THEN
335            inumov(:,:,:,:,:)=0
336            zoemea(:,:,:,:,:)=0
337            zovmea(:,:,:,:,:)=0
338         ENDIF
339         IF (nbgerr>0) THEN
340            inumbv(:,:,:,:,:)=0
341            zbemea(:,:,:,:,:)=0
342            zbvmea(:,:,:,:,:)=0
343         ENDIF
344         inuma(:,:,:,:)=0
345         zoamean(:,:,:,:)=0.0
346      ELSE
347         IF (fbdata%nvar/=nvar) THEN
348            WRITE(*,*)'Different number of nvar ',fbdata%nvar,' in ',&
349               & TRIM(filename)
350            CALL abort
351         ENDIF
352         jadd = 0
353         DO ja= 1, fbdata%nadd
354            IF (fbdata%caddname(ja)(1:2)=='Hx') THEN
355               jadd=jadd+1
356            ENDIF
357         END DO
358         IF (jadd/=nadd) THEN
359            WRITE(*,*)'Different number of nadd ',jadd,' in ',&
360               & TRIM(filename)
361            CALL abort
362         ENDIF
363         jadd = 0
364         DO ja= 1, fbdata%nadd
365            IF (fbdata%caddname(ja)(1:3)=='OBE') THEN
366               jadd=jadd+1
367            ENDIF
368         END DO
369         IF (jadd/=noberr) THEN
370            WRITE(*,*)'Different number of noberr ',jadd,' in ',&
371               & TRIM(filename)
372            CALL abort
373         ENDIF
374         jadd = 0
375         DO ja= 1, fbdata%nadd
376            IF (fbdata%caddname(ja)(1:3)=='BGE') THEN
377               jadd=jadd+1
378            ENDIF
379         END DO
380         IF (jadd/=nbgerr) THEN
381            WRITE(*,*)'Different number of nbgerr ',jadd,' in ',&
382               & TRIM(filename)
383            CALL abort
384         ENDIF
385         IF (ipdcst>0) THEN
386            IF (ipdcst>fbdata%next) THEN
387               WRITE(*,*)'Distrance to coast in file not compatible with first file'
388               CALL abort
389            ENDIF
390            IF (TRIM(fbdata%cextname(ipdcst))/='DCST') THEN
391               WRITE(*,*)'Distrance to coast in file not compatible with first file'
392               CALL abort
393            ENDIF
394         ENDIF
395      ENDIF
396      IF (lrmmean) THEN
397         CALL fb_rmmean(fbdata)
398      ENDIF
399      CALL fb_stat(fbdata,lskipbox,nmlev,zlev,lnear,nqc,nqco,&
400         &         lhistogram,lxyplot,ntyp,ctyp,ipdcst,mindcst)
401      CALL dealloc_obfbdata(fbdata)
402   ENDDO
403
404   DO jt=1,ntyp
405      DO jvar=1,nvar
406         DO jadd=1,nadd
407            jboxl=0
408            DO jbox=1,nboxuser
409               IF (lskipbox(jbox)) CYCLE
410               jboxl=jboxl+1
411               DO jlev = 1, nmlev
412                  IF ( inum(jlev,jboxl,jadd,jvar,jt) > 0 ) THEN
413                     zbias(jlev,jboxl,jadd,jvar,jt) = &
414                        & zbias(jlev,jboxl,jadd,jvar,jt)/inum(jlev,jboxl,jadd,jvar,jt)
415                     zrms(jlev,jboxl,jadd,jvar,jt) = &
416                        & SQRT(zrms(jlev,jboxl,jadd,jvar,jt)/inum(jlev,jboxl,jadd,jvar,jt))
417                     zsdev(jlev,jboxl,jadd,jvar,jt) = &
418                        & SQRT(MAX(zrms(jlev,jboxl,jadd,jvar,jt)**2-zbias(jlev,jboxl,jadd,jvar,jt)**2,0.0))
419                     zomean(jlev,jboxl,jadd,jvar,jt) = &
420                        & zomean(jlev,jboxl,jadd,jvar,jt)/inum(jlev,jboxl,jadd,jvar,jt)
421                     zmmean(jlev,jboxl,jadd,jvar,jt) = &
422                        & zmmean(jlev,jboxl,jadd,jvar,jt)/inum(jlev,jboxl,jadd,jvar,jt)
423                  ELSE
424                     zbias(jlev,jboxl,jadd,jvar,jt) = fbrmdi
425                     zrms(jlev,jboxl,jadd,jvar,jt) = fbrmdi
426                     zsdev(jlev,jboxl,jadd,jvar,jt) = fbrmdi
427                     zomean(jlev,jboxl,jadd,jvar,jt) = fbrmdi
428                     zmmean(jlev,jboxl,jadd,jvar,jt) = fbrmdi
429                  ENDIF
430               ENDDO
431            ENDDO
432         ENDDO
433         DO jadd=1,noberr
434            jboxl=0
435            DO jbox=1,nboxuser
436               IF (lskipbox(jbox)) CYCLE
437               jboxl=jboxl+1
438               DO jlev = 1, nmlev
439                  IF ( inumov(jlev,jboxl,jadd,jvar,jt) > 0 ) THEN
440                     zoemea(jlev,jboxl,jadd,jvar,jt) = &
441                        & zoemea(jlev,jboxl,jadd,jvar,jt)/inumov(jlev,jboxl,jadd,jvar,jt)
442                     zovmea(jlev,jboxl,jadd,jvar,jt) = &
443                        & zovmea(jlev,jboxl,jadd,jvar,jt)/inumov(jlev,jboxl,jadd,jvar,jt)
444                  ELSE
445                     zoemea(jlev,jboxl,jadd,jvar,jt) = fbrmdi
446                     zovmea(jlev,jboxl,jadd,jvar,jt) = fbrmdi
447                  ENDIF
448               ENDDO
449            ENDDO
450         ENDDO
451         DO jadd=1,nbgerr
452            jboxl=0
453            DO jbox=1,nboxuser
454               IF (lskipbox(jbox)) CYCLE
455               jboxl=jboxl+1
456               DO jlev = 1, nmlev
457                  IF ( inumbv(jlev,jboxl,jadd,jvar,jt) > 0 ) THEN
458                     zbemea(jlev,jboxl,jadd,jvar,jt) = &
459                        & zbemea(jlev,jboxl,jadd,jvar,jt)/inumbv(jlev,jboxl,jadd,jvar,jt)
460                     zbvmea(jlev,jboxl,jadd,jvar,jt) = &
461                        & zbvmea(jlev,jboxl,jadd,jvar,jt)/inumbv(jlev,jboxl,jadd,jvar,jt)
462                  ELSE
463                     zbemea(jlev,jboxl,jadd,jvar,jt) = fbrmdi
464                     zbvmea(jlev,jboxl,jadd,jvar,jt) = fbrmdi
465                  ENDIF
466               ENDDO
467            ENDDO
468         ENDDO
469      ENDDO
470   ENDDO
471   DO jt=1,ntyp
472      DO jvar=1,nvar
473         jboxl=0
474         DO jbox=1,nboxuser
475            IF (lskipbox(jbox)) CYCLE
476            jboxl=jboxl+1
477            DO jlev = 1, nmlev
478               IF ( inuma(jlev,jboxl,jvar,jt) > 0 ) THEN
479                  zoamean(jlev,jboxl,jvar,jt) = &
480                     & zoamean(jlev,jboxl,jvar,jt)/inuma(jlev,jboxl,jvar,jt)
481               ELSE
482                  zoamean(jlev,jboxl,jvar,jt) = fbrmdi
483               ENDIF
484            ENDDO
485         ENDDO
486      ENDDO
487   ENDDO
488
489   IF (ltext) THEN
490
491      DO jt=1,ntyp
492         DO jvar=1,nvar
493            DO jadd=1,nadd
494               jboxl=0
495               DO jbox=1,nboxuser
496                  IF (lskipbox(jbox)) CYCLE
497                  jboxl=jboxl+1
498                  WRITE(filename,'(7A)')TRIM(cname(jvar)),&
499                     &                  TRIM(caddname(jadd)),'_',&
500                     &                  TRIM(cl_boxes_user(jbox)),'_',&
501                     &                  TRIM(ADJUSTL(ctyp(jt))),'.dat'
502                  OPEN(10,file=TRIM(filename))
503                  DO jlev = 1, nmlev
504                     WRITE(10,'(F16.7,2I12,5F17.10)') zlev(jlev), &
505                        & jlev, inum(jlev,jboxl,jadd,jvar,jt), &
506                        & zbias(jlev,jboxl,jadd,jvar,jt), &
507                        & zrms(jlev,jboxl,jadd,jvar,jt), &
508                        & zsdev(jlev,jboxl,jadd,jvar,jt), &
509                        & zomean(jlev,jboxl,jadd,jvar,jt), &
510                        & zmmean(jlev,jboxl,jadd,jvar,jt)
511                  ENDDO
512                  CLOSE(10)
513               ENDDO
514            ENDDO
515            DO jadd=1,noberr
516               jboxl=0
517               DO jbox=1,nboxuser
518                  IF (lskipbox(jbox)) CYCLE
519                  jboxl=jboxl+1
520                  WRITE(filename,'(7A)')TRIM(cname(jvar)),&
521                     &                  TRIM(cobename(jadd)),'_',&
522                     &                  TRIM(cl_boxes_user(jbox)),'_',&
523                     &                  TRIM(ADJUSTL(ctyp(jt))),'.dat'
524                  OPEN(10,file=TRIM(filename))
525                  DO jlev = 1, nmlev
526                     WRITE(10,'(F16.7,2I12,5F17.10)') zlev(jlev), &
527                        & jlev, inumov(jlev,jboxl,jadd,jvar,jt), &
528                        & zoemea(jlev,jboxl,jadd,jvar,jt), &
529                        & zovmea(jlev,jboxl,jadd,jvar,jt)
530                  ENDDO
531                  CLOSE(10)
532               ENDDO
533            ENDDO
534            DO jadd=1,nbgerr
535               jboxl=0
536               DO jbox=1,nboxuser
537                  IF (lskipbox(jbox)) CYCLE
538                  jboxl=jboxl+1
539                  WRITE(filename,'(7A)')TRIM(cname(jvar)),&
540                     &                  TRIM(cbgename(jadd)),'_',&
541                     &                  TRIM(cl_boxes_user(jbox)),'_',&
542                     &                  TRIM(ADJUSTL(ctyp(jt))),'.dat'
543                  OPEN(10,file=TRIM(filename))
544                  DO jlev = 1, nmlev
545                     WRITE(10,'(F16.7,2I12,5F17.10)') zlev(jlev), &
546                        & jlev, inumbv(jlev,jboxl,jadd,jvar,jt), &
547                        & zbemea(jlev,jboxl,jadd,jvar,jt), &
548                        & zbvmea(jlev,jboxl,jadd,jvar,jt)
549                  ENDDO
550                  CLOSE(10)
551               ENDDO
552            ENDDO
553         ENDDO
554      ENDDO
555
556      DO jt=1,ntyp
557         DO jvar=1,nvar
558            jboxl=0
559            DO jbox=1,nboxuser
560               IF (lskipbox(jbox)) CYCLE
561               jboxl=jboxl+1
562               WRITE(filename,'(7A)')TRIM(cname(jvar)),'_',&
563                  &                  TRIM(cl_boxes_user(jbox)),'_',&
564                  &                  TRIM(ADJUSTL(ctyp(jt))),'.dat'
565               OPEN(10,file=TRIM(filename))
566               DO jlev = 1, nmlev
567                  WRITE(10,'(F16.7,2I12,F17.10)') zlev(jlev), &
568                     & jlev, inuma(jlev,jboxl,jvar,jt), &
569                     & zoamean(jlev,jboxl,jvar,jt)
570               ENDDO
571               CLOSE(10)
572            ENDDO
573         ENDDO
574      ENDDO
575
576      IF (lhistogram) THEN
577         DO jt=1,ntyp
578            DO jvar=1,nvar
579               DO jadd=1,nadd
580                  jboxl=0
581                  DO jbox=1,nboxuser
582                     IF (lskipbox(jbox)) CYCLE
583                     jboxl=jboxl+1
584                     WRITE(filename,'(7A)')TRIM(cname(jvar)),&
585                        &                  TRIM(caddname(jadd)),'_',&
586                        &                  TRIM(cl_boxes_user(jbox)),'_',&
587                        &                  TRIM(ADJUSTL(ctyp(jt))),&
588                        &                  '_histogram.dat'
589                     OPEN(10,file=TRIM(filename))
590                     WRITE(10,'(A10,1000F10.2)')'#',(zlev(jlev),jlev=1,nmlev)
591                     DO ji=1,hist(jvar)%npoints
592                        WRITE(10,'(F10.2,1000I10)') &
593                           & zhistmin(jvar)+(ji-1)*zhiststep(jvar), &
594                           & (hist(jvar)%nhist(ji,jlev,jboxl,jadd,jt),jlev=1,nmlev)
595                     ENDDO
596                     CLOSE(10)
597                  ENDDO
598               ENDDO
599            ENDDO
600         ENDDO
601      ENDIF
602     
603   ENDIF
604
605   IF (lnetcdf) THEN
606
607      IF (nadd>0) THEN
608         DO jt=1,ntyp
609            WRITE(outfilename,'(3A)')'fbstat_',TRIM(ADJUSTL(ctyp(jt))),'.nc'
610            CALL fbstat_ncwrite(TRIM(outfilename),&
611               & nvar,cname,nadd,caddname,noberr,cobename,nbgerr,cbgename,&
612               & nboxuser,nboxuserl,20,cl_boxes_user,lskipbox,nmlev,zlev,&
613               & inum(:,:,:,:,jt),zbias(:,:,:,:,jt),zrms(:,:,:,:,jt), &
614               & zsdev(:,:,:,:,jt),zomean(:,:,:,:,jt),zmmean(:,:,:,:,jt),&
615               & inuma(:,:,:,jt),zoamean(:,:,:,jt), &
616               & inumov(:,:,:,:,jt),zoemea(:,:,:,:,jt),zovmea(:,:,:,:,jt), &
617               & inumbv(:,:,:,:,jt),zbemea(:,:,:,:,jt),zbvmea(:,:,:,:,jt) )
618            IF (lhistogram) THEN
619               WRITE(outfilename,'(3A)')'fbstat_hist_',TRIM(ADJUSTL(ctyp(jt))),'.nc'
620               CALL fbstat_ncwrite_hist(TRIM(outfilename),&
621                  & nvar,cname,nadd,caddname,&
622                  & nboxuser,20,cl_boxes_user,lskipbox,nmlev,zlev,&
623                  & hist,zhistmin,zhiststep,jt)
624            ENDIF
625            IF (lxyplot) THEN
626               WRITE(outfilename,'(3A)')'fbstat_xyplot_',TRIM(ADJUSTL(ctyp(jt))),'.nc'
627               CALL fbstat_ncwrite_xy(TRIM(outfilename),&
628                  & nvar,cname,nadd,caddname,&
629                  & nboxuser,20,cl_boxes_user,lskipbox,nmlev,zlev,&
630                  & xy,zxymin,zxystep,jt)
631            ENDIF
632         ENDDO
633      ENDIF
634   ENDIF
635
636   IF (lomona) THEN
637
638      IF (ntyp>1) THEN
639         WRITE(*,*)'Omona file only supported for the first type which is : ',TRIM(ctyp(1))
640      ENDIF
641      IF (nmlev>1) THEN
642         ALLOCATE(zdat3d(nmlev,nboxuser,1))
643      ELSE
644         ALLOCATE(zdat2d(nboxuser,1))
645      ENDIF
646
647      cl_expnam=expver
648      WRITE(cl_date,'(I8.8)')inidate
649      i_dp = nmlev
650      itime=icurdate
651      linnerp=.TRUE.
652      iloopno = loopno
653      linnerini = linner
654      i_fill=0
655
656      DO jt=1,ntyp
657         DO jvar = 1, nvar
658            linner = linnerini
659            loopno = iloopno
660            SELECT CASE (TRIM(cname(jvar)))
661            CASE('POTM')
662               cl_var = 'votemper'
663            CASE('PSAL')
664               cl_var='vosaline'
665            CASE('SLA')
666               cl_var='sossheig'
667            CASE('SST')
668               cl_var='sosstsst'
669            CASE('UVEL')
670               cl_var='vozocrtx'
671            CASE('VVEL')
672               cl_var='vomecrty'
673            END SELECT
674            DO jadd = 1, nadd
675               linner = (caddname(jadd)(1:3)=='Hxa').OR.linner
676               IF (lpassive) THEN
677                  ityp=145
678               ELSE
679                  IF (linner) THEN
680                     linnerp=.TRUE.
681                     ityp=144
682                     IF (jadd>1) loopno=loopno+1
683                  ELSE
684                     ityp=142
685                     IF (.NOT.linnerp) THEN
686                        IF (jadd>1) loopno=loopno+1
687                     ENDIF
688                  ENDIF
689               ENDIF
690               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
691                  & TRIM(ADJUSTL(ctyp(jt)))
692               CALL obs_variable_att(cltyp)
693               IF (nmlev>1) THEN
694                  zdat3d(:,:,1) = zbias(:,:,jadd,jvar,jt)
695                  i_fill=0
696                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
697                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
698                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
699               ELSE
700                  zdat2d(:,1) = zbias(1,:,jadd,jvar,jt)
701                  i_fill=0
702                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
703                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
704               ENDIF
705               IF (lpassive) THEN
706                  ityp=245
707               ELSE
708                  IF (linner) THEN
709                     ityp=244
710                  ELSE
711                     ityp=242
712                  ENDIF
713               ENDIF
714               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
715                  & TRIM(ADJUSTL(ctyp(jt)))
716               CALL obs_variable_att(cltyp)
717               IF (nmlev>1) THEN
718                  zdat3d(:,:,1) = zrms(:,:,jadd,jvar,jt)
719                  i_fill=0
720                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
721                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
722                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
723               ELSE
724                  zdat2d(:,1) = zrms(1,:,jadd,jvar,jt)
725                  i_fill=0
726                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
727                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
728               ENDIF
729               IF (lpassive) THEN
730                  ityp=345
731               ELSE
732                  IF (linner) THEN
733                     ityp=344
734                  ELSE
735                     ityp=342
736                  ENDIF
737               ENDIF
738               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
739               & TRIM(ADJUSTL(ctyp(jt)))
740               CALL obs_variable_att(cltyp)
741               IF (nmlev>1) THEN
742                  zdat3d(:,:,1) = zsdev(:,:,jadd,jvar,jt)
743                  i_fill=0
744                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
745                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
746                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
747               ELSE
748                  zdat2d(:,1) = zsdev(1,:,jadd,jvar,jt)
749                  i_fill=0
750                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
751                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
752               ENDIF
753               IF (lpassive) THEN
754                  ityp=445
755               ELSE
756                  IF (linner) THEN
757                     ityp=444
758                  ELSE
759                     ityp=442
760                  ENDIF
761               ENDIF
762               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
763                  & TRIM(ADJUSTL(ctyp(jt)))
764               CALL obs_variable_att(cltyp)
765               IF (nmlev>1) THEN
766                  zdat3d(:,:,1) = inum(:,:,jadd,jvar,jt)
767                  i_fill=0
768                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
769                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
770                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
771               ELSE
772                  zdat2d(:,1) = inum(1,:,jadd,jvar,jt)
773                  i_fill=0
774                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
775                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
776               ENDIF
777               IF (lpassive) THEN
778                  ityp=545
779               ELSE
780                  IF (linner) THEN
781                     ityp=544
782                  ELSE
783                     ityp=542
784                  ENDIF
785               ENDIF
786               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
787                  & TRIM(ADJUSTL(ctyp(jt)))
788               CALL obs_variable_att(cltyp)
789               IF (nmlev>1) THEN
790                  zdat3d(:,:,1) = zomean(:,:,jadd,jvar,jt)
791                  i_fill=0
792                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
793                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
794                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
795               ELSE
796                  zdat2d(:,1) = zomean(1,:,jadd,jvar,jt)
797                  i_fill=0
798                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
799                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
800               ENDIF
801               IF (lpassive) THEN
802                  ityp=645
803               ELSE
804                  IF (linner) THEN
805                     ityp=644
806                  ELSE
807                     ityp=642
808                  ENDIF
809               ENDIF
810               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
811                  & TRIM(ADJUSTL(ctyp(jt)))
812               CALL obs_variable_att(cltyp)
813               IF (nmlev>1) THEN
814                  zdat3d(:,:,1) = zmmean(:,:,jadd,jvar,jt)
815                  i_fill=0
816                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
817                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
818                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
819               ELSE
820                  zdat2d(:,1) = zmmean(1,:,jadd,jvar,jt)
821                  i_fill=0
822                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
823                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
824               ENDIF
825               linner=.FALSE.
826            ENDDO
827            loopno = iloopno
828            DO jadd = 1, noberr
829               linner = .TRUE.
830               ityp = 139
831               IF (jadd>1) loopno=loopno+1
832               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
833                  & TRIM(ADJUSTL(ctyp(jt)))
834               CALL obs_variable_att(cltyp)
835               IF (nmlev>1) THEN
836                  zdat3d(:,:,1) = zoemea(:,:,jadd,jvar,jt)
837                  i_fill=0
838                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
839                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
840                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
841               ELSE
842                  zdat2d(:,1) = zoemea(1,:,jadd,jvar,jt)
843                  i_fill=0
844                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
845                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
846               ENDIF
847               ityp = 239
848               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
849                  & TRIM(ADJUSTL(ctyp(jt)))
850               CALL obs_variable_att(cltyp)
851               IF (nmlev>1) THEN
852                  zdat3d(:,:,1) = zovmea(:,:,jadd,jvar,jt)
853                  i_fill=0
854                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
855                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
856                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
857               ELSE
858                  zdat2d(:,1) = zovmea(1,:,jadd,jvar,jt)
859                  i_fill=0
860                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
861                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
862               ENDIF
863            ENDDO
864            loopno = iloopno
865            DO jadd = 1, nbgerr
866               linner = .TRUE.
867               ityp = 141
868               IF (jadd>1) loopno=loopno+1
869               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
870                  & TRIM(ADJUSTL(ctyp(jt)))
871               CALL obs_variable_att(cltyp)
872               IF (nmlev>1) THEN
873                  zdat3d(:,:,1) = zbemea(:,:,jadd,jvar,jt)
874                  i_fill=0
875                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
876                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
877                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
878               ELSE
879                  zdat2d(:,1) = zbemea(1,:,jadd,jvar,jt)
880                  i_fill=0
881                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
882                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
883               ENDIF
884               ityp = 241
885               WRITE(cltyp,'(I3.3,A1,I2.2,A1,A)')ityp,'_',loopno,'_',&
886                  & TRIM(ADJUSTL(ctyp(jt)))
887               CALL obs_variable_att(cltyp)
888               IF (nmlev>1) THEN
889                  zdat3d(:,:,1) = zbvmea(:,:,jadd,jvar,jt)
890                  i_fill=0
891                  CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
892                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
893                  CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
894               ELSE
895                  zdat2d(:,1) = zbvmea(1,:,jadd,jvar,jt)
896                  i_fill=0
897                  CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
898                     &                    cl_boxes_user,REAL(fbrmdi),i_fill)
899               ENDIF
900            ENDDO
901            IF (lpassive) THEN
902               ityp=745
903            ELSE
904               ityp=742
905            ENDIF
906            WRITE(cltyp,'(I3.3,A1,A)')ityp,'_',&
907               & TRIM(ADJUSTL(ctyp(jt)))
908            CALL obs_variable_att(cltyp)
909            IF (nmlev>1) THEN
910               zdat3d(:,:,1) = inuma(:,:,jvar,jt)
911               i_fill=0
912               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
913                  &                    cl_boxes_user,REAL(fbrmdi),i_fill)
914               CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
915            ELSE
916               zdat2d(:,1) = inuma(1,:,jvar,jt)
917               i_fill=0
918               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
919                  &                    cl_boxes_user,REAL(fbrmdi),i_fill)
920            ENDIF
921            IF (lpassive) THEN
922               ityp=845
923            ELSE
924               ityp=842
925            ENDIF
926            WRITE(cltyp,'(I3.3,A1,A)')ityp,'_',&
927               & TRIM(ADJUSTL(ctyp(jt)))
928            CALL obs_variable_att(cltyp)
929            IF (nmlev>1) THEN
930               zdat3d(:,:,1) = zoamean(:,:,jvar,jt)
931               i_fill=0
932               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
933                  &                    cl_boxes_user,REAL(fbrmdi),i_fill)
934               CALL write_dep_netcdf(cl_filename_out,cl_boxes_user,zlev)
935            ELSE
936               zdat2d(:,1) = zoamean(1,:,jvar,jt)
937               i_fill=0
938               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
939                  &                    cl_boxes_user,REAL(fbrmdi),i_fill)
940            ENDIF
941         ENDDO
942      ENDDO
943
944      IF (nmlev>1) THEN
945         DEALLOCATE(zdat3d)
946      ELSE
947         DEALLOCATE(zdat2d)
948      ENDIF
949     
950   ENDIF
951
952CONTAINS
953
954   SUBROUTINE fb_stat(fbdata,lskipbox,nmlev,zlev,lnear,kqc,kqco,lhist,lxyplot,&
955      &               ntyp,ctyp,ipdcst,mindcst)
956      USE fbaccdata
957      USE coords
958      TYPE(obfbdata) :: fbdata
959      LOGICAL, DIMENSION(nboxuser) :: lskipbox
960      INTEGER :: nmlev
961      REAL :: zlev(nmlev)
962      LOGICAL :: lnear
963      INTEGER :: kqc,kqco
964      LOGICAL :: lhist,lxyplot
965      INTEGER :: ntyp
966      CHARACTER(len=ilentyp), DIMENSION(ntyp) :: ctyp
967      INTEGER :: ipdcst
968      REAL :: mindcst
969      INTEGER, DIMENSION(nboxuser) :: jlboxnum
970      INTEGER :: jlev, jobs, jvar, klev,jlev2,ih,ja,jadd,jbox,jt,ix,iy,jboxl
971      REAL :: zarea(4),zlat,zlon,zdiff,zdiff2,zvar
972     
973      jboxl=0
974      jlboxnum=-1
975      DO jbox = 1, nboxuser
976         IF (lskipbox(jbox)) CYCLE
977         jboxl=jboxl+1
978         jlboxnum(jbox)=jboxl
979      ENDDO
980
981      !$omp parallel do default(shared) private(jlev,jobs,jvar,klev,jlev2,ih,ja,jadd,jbox,jt,ix,iy,jboxl,zarea,zlat,zlon,zdiff,zdiff2)
982      DO jbox = 1, nboxuser
983         IF (lskipbox(jbox)) CYCLE
984         jboxl=jlboxnum(jbox)
985         CALL coord_area_user(cl_boxes_user(jbox),zarea)
986         DO jobs = 1, fbdata%nobs
987            ! Reject observations with observation, position or time flag rejections
988            IF (fbdata%ioqc(jobs)>kqco) CYCLE
989            IF (fbdata%ipqc(jobs)>kqco) CYCLE
990            IF (fbdata%itqc(jobs)>kqco) CYCLE
991            zlat = fbdata%pphi(jobs)
992            zlon = fbdata%plam(jobs)
993            IF (zlon<0) zlon=zlon+360
994            IF (zlon>360) zlon=zlon-360
995            IF ( ( zlat .GE. zarea(3) ) .AND. &
996               & ( zlat .LE. zarea(4) ) .AND. &
997               & ( ( ( zlon .GE. zarea(1) ) .AND. &
998               &     ( zlon .LE. zarea(2) ) ) .OR. &
999               &   ( ( zarea(2) .LE. zarea(1) ) .AND. &
1000               &     ( zlon .GE. zarea(1) ) .AND. &
1001               &     ( zlon .LE. 360        ) ) .OR. &
1002               &   ( ( zarea(2)   .LE. zarea(1) ) .AND. &
1003               &     ( zlon .GE. 0          ) .AND. &
1004               &     ( zlon .LE. zarea(2) ) ) ) ) THEN
1005           
1006               DO jlev = 1, fbdata%nlev
1007                  IF (ipdcst>0) THEN
1008                     IF (fbdata%pext(jlev,jobs,ipdcst)==fbrmdi) CYCLE
1009                     IF (fbdata%pext(jlev,jobs,ipdcst)<mindcst) CYCLE
1010                  ENDIF
1011                  DO jvar = 1, fbdata%nvar
1012                     IF (nmlev==1) THEN
1013                        klev=1
1014                     ELSE
1015                        IF (lnear) THEN
1016                           zdiff=ABS(fbdata%pdep(jlev,jobs)-zlev(1))
1017                           klev=1
1018                           DO jlev2=2,nmlev
1019                              zdiff2=ABS(fbdata%pdep(jlev,jobs)-zlev(jlev2))
1020                              IF (zdiff2<zdiff) THEN
1021                                 klev=jlev2
1022                                 zdiff=zdiff2
1023                              ENDIF
1024                           ENDDO
1025                        ELSE
1026                           klev = fbdata%iobsk(jlev,jobs,jvar)-1
1027                        ENDIF
1028                        IF ( klev > nmlev ) THEN
1029                           DO ja = 1, fbdata%nadd
1030                              IF ( fbdata%caddname(ja)(1:2) /= 'Hx' ) CYCLE
1031                              IF ( ABS(fbdata%padd(jlev,jobs,ja,jvar))<9000 ) THEN
1032                                 WRITE(*,*)'Error in fb_stat'
1033                                 WRITE(*,*)'Increase nmlev to at least ',klev
1034                                 klev=nmlev
1035                                 CALL abort
1036                              ENDIF
1037                           ENDDO
1038                        ENDIF
1039                     ENDIF
1040                     IF (( klev > 0 ).AND. &
1041                        &(ABS(fbdata%pob(jlev,jobs,jvar)) < 9000 )) THEN
1042                        DO jt=1,ntyp
1043                           IF (TRIM(ADJUSTL(ctyp(jt)))/='all') THEN
1044                              IF (TRIM(ADJUSTL(ctyp(jt)))/=TRIM(ADJUSTL(fbdata%cdtyp(jobs)))) CYCLE
1045                           ENDIF
1046                           inuma(klev,jboxl,jvar,jt) = inuma(klev,jboxl,jvar,jt) + 1
1047                           zoamean(klev,jboxl,jvar,jt) = zoamean(klev,jboxl,jvar,jt) + &
1048                              & fbdata%pob(jlev,jobs,jvar)
1049                        ENDDO
1050                     ENDIF
1051                     IF ( fbdata%ivlqc(jlev,jobs,jvar) < 0 ) CYCLE
1052                     IF ( fbdata%ivlqc(jlev,jobs,jvar) > kqc ) CYCLE
1053                     IF (( klev > 0 ).AND. &
1054                        &(ABS(fbdata%pob(jlev,jobs,jvar)) < 9000 )) THEN
1055                        jadd = 0
1056                        DO ja = 1, fbdata%nadd
1057                           IF ( fbdata%caddname(ja)(1:2) /= 'Hx' ) CYCLE
1058                           jadd = jadd + 1
1059                           IF ( ABS(fbdata%padd(jlev,jobs,ja,jvar)) < 9000 ) THEN
1060                              zdiff = ( fbdata%padd(jlev,jobs,ja,jvar) - &
1061                                 &      fbdata%pob(jlev,jobs,jvar) )
1062                              DO jt=1,ntyp
1063                                 IF (TRIM(ADJUSTL(ctyp(jt)))/='all') THEN
1064                                    IF (TRIM(ADJUSTL(ctyp(jt)))/=TRIM(ADJUSTL(fbdata%cdtyp(jobs)))) CYCLE
1065                                 ENDIF
1066                                 inum(klev,jboxl,jadd,jvar,jt) = inum(klev,jboxl,jadd,jvar,jt) + 1
1067                                 zbias(klev,jboxl,jadd,jvar,jt) = zbias(klev,jboxl,jadd,jvar,jt) + &
1068                                    & zdiff
1069                                 zrms(klev,jboxl,jadd,jvar,jt) = zrms(klev,jboxl,jadd,jvar,jt) + &
1070                                    & zdiff * zdiff
1071                                 zomean(klev,jboxl,jadd,jvar,jt) = zomean(klev,jboxl,jadd,jvar,jt) + &
1072                                    & fbdata%pob(jlev,jobs,jvar)
1073                                 zmmean(klev,jboxl,jadd,jvar,jt) = zmmean(klev,jboxl,jadd,jvar,jt) + &
1074                                    & fbdata%padd(jlev,jobs,ja,jvar)
1075                              ENDDO
1076                              IF (ABS(zdiff)>zcheck(jvar)) THEN
1077                                 WRITE(*,*)'Departure outside check range ',&
1078                                    & TRIM(fbdata%cname(jvar)),' entry ',&
1079                                    & fbdata%caddname(jadd)
1080                                 WRITE(*,*)'Depar = ',zdiff
1081                                 WRITE(*,*)'Check = ',zcheck(jvar)
1082                                 WRITE(*,*)'Id    = ',fbdata%cdwmo(jobs)
1083                                 WRITE(*,*)'Lat   = ',fbdata%pphi(jobs)
1084                                 WRITE(*,*)'Lon   = ',fbdata%plam(jobs)
1085                                 WRITE(*,*)'Tim   = ',fbdata%ptim(jobs)
1086                                 WRITE(*,*)'Depth = ',fbdata%pdep(jlev,jobs)
1087                                 WRITE(*,*)'Obs   = ',fbdata%pob(jlev,jobs,jvar)
1088                                 WRITE(*,*)'Var   = ',fbdata%padd(jlev,jobs,ja,jvar)
1089                                 WRITE(*,*)'QC    = ',fbdata%ivlqc(jlev,jobs,jvar)
1090                                 WRITE(*,*)'QCflag= ',fbdata%ivlqcf(:,jlev,jobs,jvar)
1091                              ENDIF
1092                              IF (lhist) THEN
1093                                 ih=NINT((zdiff-zhistmin(jvar))/zhiststep(jvar))+1
1094                                 IF ((ih>=1).AND.(ih<=hist(jvar)%npoints)) THEN
1095                                    DO jt=1,ntyp
1096                                       IF (TRIM(ADJUSTL(ctyp(jt)))/='all') THEN
1097                                          IF (TRIM(ADJUSTL(ctyp(jt)))/=TRIM(ADJUSTL(fbdata%cdtyp(jobs)))) CYCLE
1098                                       ENDIF
1099                                       hist(jvar)%nhist(ih,klev,jboxl,jadd,jt) = &
1100                                          hist(jvar)%nhist(ih,klev,jboxl,jadd,jt) +1
1101                                    ENDDO
1102                                 ELSE
1103                                    WRITE(*,*)'Histogram value outside range for ',&
1104                                       & TRIM(fbdata%cname(jvar)),' entry ',&
1105                                       & fbdata%caddname(jadd)
1106                                    WRITE(*,*)'Value = ',zdiff
1107                                    WRITE(*,*)'Range = ',zhistmin(jvar),zhistmax(jvar)
1108                                    WRITE(*,*)'Step  = ',zhiststep(jvar)
1109                                    WRITE(*,*)'Index = ',ih
1110                                    WRITE(*,*)'Id    = ',TRIM(fbdata%cdwmo(jobs))
1111                                    WRITE(*,*)'Typ   = ',TRIM(fbdata%cdtyp(jobs))
1112                                    WRITE(*,*)'Lat   = ',fbdata%pphi(jobs)
1113                                    WRITE(*,*)'Lon   = ',fbdata%plam(jobs)
1114                                    WRITE(*,*)'Tim   = ',fbdata%ptim(jobs)
1115                                    WRITE(*,*)'Depth = ',fbdata%pdep(jlev,jobs)
1116                                    WRITE(*,*)'Obs   = ',fbdata%pob(jlev,jobs,jvar)
1117                                    WRITE(*,*)'Var   = ',fbdata%padd(jlev,jobs,jadd,jvar)
1118                                    WRITE(*,*)'QC    = ',fbdata%ivlqc(jlev,jobs,jvar)
1119                                    WRITE(*,*)'QCflag= ',fbdata%ivlqcf(:,jlev,jobs,jvar)
1120                                 ENDIF
1121                              ENDIF
1122                              IF (lxyplot) THEN
1123                                 ix=NINT((fbdata%pob(jlev,jobs,jvar)-zxymin(jvar))/&
1124                                    &    zxystep(jvar))+1
1125                                 iy=NINT((fbdata%padd(jlev,jobs,ja,jvar)-zxymin(jvar))/&
1126                                    &    zxystep(jvar))+1
1127                                 IF ((ix>=1).AND.(ix<=xy(jvar)%npoints).AND. &
1128                                    &(iy>=1).AND.(iy<=xy(jvar)%npoints)) THEN
1129                                    DO jt=1,ntyp
1130                                       IF (TRIM(ADJUSTL(ctyp(jt)))/='all') THEN
1131                                          IF (TRIM(ADJUSTL(ctyp(jt)))/=TRIM(ADJUSTL(fbdata%cdtyp(jobs)))) CYCLE
1132                                       ENDIF
1133                                       xy(jvar)%nxy(ix,iy,klev,jboxl,jadd,jt) = &
1134                                          xy(jvar)%nxy(ix,iy,klev,jboxl,jadd,jt) +1
1135                                    ENDDO
1136                                 ELSE
1137                                    WRITE(*,*)'xy plot values outside range for ',&
1138                                       & TRIM(fbdata%cname(jvar)),' entry ',&
1139                                       & fbdata%caddname(jadd)
1140                                    WRITE(*,*)'Obs   = ',fbdata%pob(jlev,jobs,jvar)
1141                                    WRITE(*,*)'Model = ',fbdata%padd(jlev,jobs,ja,jvar)
1142                                    WRITE(*,*)'Range = ',zxymin(jvar),zxymax(jvar)
1143                                    WRITE(*,*)'Step  = ',zxystep(jvar)
1144                                    WRITE(*,*)'Index = ',ih
1145                                    WRITE(*,*)'Id    = ',TRIM(fbdata%cdwmo(jobs))
1146                                    WRITE(*,*)'Typ   = ',TRIM(fbdata%cdtyp(jobs))
1147                                    WRITE(*,*)'Lat   = ',fbdata%pphi(jobs)
1148                                    WRITE(*,*)'Lon   = ',fbdata%plam(jobs)
1149                                    WRITE(*,*)'Tim   = ',fbdata%ptim(jobs)
1150                                    WRITE(*,*)'Depth = ',fbdata%pdep(jlev,jobs)
1151                                    WRITE(*,*)'Obs   = ',fbdata%pob(jlev,jobs,jvar)
1152                                    WRITE(*,*)'Var   = ',fbdata%padd(jlev,jobs,jadd,jvar)
1153                                    WRITE(*,*)'QC    = ',fbdata%ivlqc(jlev,jobs,jvar)
1154                                    WRITE(*,*)'QCflag= ',fbdata%ivlqcf(:,jlev,jobs,jvar)
1155                                 ENDIF
1156                              ENDIF
1157                           ENDIF
1158                        ENDDO
1159                        jadd = 0
1160                        DO ja = 1, fbdata%nadd
1161                           IF ( fbdata%caddname(ja)(1:3) /= 'OBE' ) CYCLE
1162                           jadd = jadd + 1
1163                           IF ( ABS(fbdata%padd(jlev,jobs,ja,jvar)) < 9000 ) THEN
1164                              zvar = fbdata%padd(jlev,jobs,ja,jvar)*fbdata%padd(jlev,jobs,ja,jvar)
1165                              DO jt=1,ntyp
1166                                 IF (TRIM(ADJUSTL(ctyp(jt)))/='all') THEN
1167                                    IF (TRIM(ADJUSTL(ctyp(jt)))/=TRIM(ADJUSTL(fbdata%cdtyp(jobs)))) CYCLE
1168                                 ENDIF
1169                                 inumov(klev,jboxl,jadd,jvar,jt) = inumov(klev,jboxl,jadd,jvar,jt) + 1
1170                                 zoemea(klev,jboxl,jadd,jvar,jt) = zoemea(klev,jboxl,jadd,jvar,jt) + &
1171                                    & fbdata%padd(jlev,jobs,ja,jvar)                                 
1172                                 zovmea(klev,jboxl,jadd,jvar,jt) = zovmea(klev,jboxl,jadd,jvar,jt) + zvar
1173                              ENDDO
1174                           ENDIF
1175                        ENDDO
1176                        jadd = 0
1177                        DO ja = 1, fbdata%nadd
1178                           IF ( fbdata%caddname(ja)(1:3) /= 'BGE' ) CYCLE
1179                           jadd = jadd + 1
1180                           IF ( ABS(fbdata%padd(jlev,jobs,ja,jvar)) < 9000 ) THEN
1181                              zvar = fbdata%padd(jlev,jobs,ja,jvar)*fbdata%padd(jlev,jobs,ja,jvar)
1182                              DO jt=1,ntyp
1183                                 IF (TRIM(ADJUSTL(ctyp(jt)))/='all') THEN
1184                                    IF (TRIM(ADJUSTL(ctyp(jt)))/=TRIM(ADJUSTL(fbdata%cdtyp(jobs)))) CYCLE
1185                                 ENDIF
1186                                 inumbv(klev,jboxl,jadd,jvar,jt) = inumbv(klev,jboxl,jadd,jvar,jt) + 1
1187                                 zbemea(klev,jboxl,jadd,jvar,jt) = zbemea(klev,jboxl,jadd,jvar,jt) + &
1188                                    & fbdata%padd(jlev,jobs,ja,jvar)
1189                                 zbvmea(klev,jboxl,jadd,jvar,jt) = zbvmea(klev,jboxl,jadd,jvar,jt) + zvar
1190                              ENDDO
1191                           ENDIF
1192                        ENDDO
1193                     ENDIF
1194                  ENDDO
1195               ENDDO
1196            ENDIF
1197         ENDDO
1198      ENDDO
1199      !$omp end parallel do
1200     
1201   END SUBROUTINE fb_stat
1202
1203   SUBROUTINE fb_rmmean(fbdata)
1204      TYPE(obfbdata) :: fbdata
1205      INTEGER :: jadd,jmean
1206
1207      jmean=0
1208      DO jadd=1,fbdata%nadd
1209         IF (fbdata%caddname(jadd)(1:4)=='MEAN') THEN
1210            jmean=jadd
1211            EXIT
1212         ENDIF
1213      ENDDO
1214      IF (jmean==0) THEN
1215         WRITE(*,*)'Warning: MEAN additional variable not found'
1216         RETURN
1217      ENDIF
1218      IF (fbdata%nobs>0) THEN
1219         DO jadd=1,fbdata%nadd 
1220            IF (fbdata%caddname(jadd)(1:2)=='Hx') THEN
1221               fbdata%padd(:,:,jadd,:)=fbdata%padd(:,:,jadd,:)&
1222                  &                   +fbdata%padd(:,:,jmean,:)
1223            ENDIF
1224         ENDDO
1225      ENDIF
1226
1227   END SUBROUTINE fb_rmmean
1228
1229   SUBROUTINE getlevsmean(nmlev,zlev)
1230      IMPLICIT NONE
1231      INTEGER :: nmlev
1232      REAL,DIMENSION(nmlev) :: zlev
1233      REAL,DIMENSION(nmlev+1) :: ztmp
1234      INTEGER :: i
1235
1236      zlev(:)=9999.9
1237      CALL getlevs(nmlev+1,ztmp)
1238      DO i=1,nmlev
1239         zlev(i)=0.5*(ztmp(i)+ztmp(i+1))
1240      ENDDO
1241
1242   END SUBROUTINE getlevsmean
1243
1244   SUBROUTINE getlevs(nmlev,zlev)
1245      IMPLICIT NONE
1246      INTEGER :: nmlev
1247      REAL,DIMENSION(nmlev) :: zlev
1248     
1249      zlev(:)=9999.9
1250
1251      IF (nmlev==42) THEN
1252         zlev(1)=5.02159
1253         zlev(2)=15.07854
1254         zlev(3)=25.16046
1255         zlev(4)=35.27829
1256         zlev(5)=45.44776
1257         zlev(6)=55.69149
1258         zlev(7)=66.04198
1259         zlev(8)=76.54591
1260         zlev(9)=87.27029
1261         zlev(10)=98.31118
1262         zlev(11)=109.8062
1263         zlev(12)=121.9519
1264         zlev(13)=135.0285
1265         zlev(14)=149.4337
1266         zlev(15)=165.7285
1267         zlev(16)=184.6975
1268         zlev(17)=207.4254
1269         zlev(18)=235.3862
1270         zlev(19)=270.5341
1271         zlev(20)=315.3741
1272         zlev(21)=372.9655
1273         zlev(22)=446.8009
1274         zlev(23)=540.5022
1275         zlev(24)=657.3229
1276         zlev(25)=799.5496
1277         zlev(26)=967.9958
1278         zlev(27)=1161.806
1279         zlev(28)=1378.661
1280         zlev(29)=1615.291
1281         zlev(30)=1868.071
1282         zlev(31)=2133.517
1283         zlev(32)=2408.583
1284         zlev(33)=2690.780
1285         zlev(34)=2978.166
1286         zlev(35)=3269.278
1287         zlev(36)=3563.041
1288         zlev(37)=3858.676
1289         zlev(38)=4155.628
1290         zlev(39)=4453.502
1291         zlev(40)=4752.021
1292         zlev(41)=5050.990
1293         zlev(42)=5350.272
1294      ELSEIF (nmlev==31) THEN
1295         zlev(1)=4.999938
1296         zlev(2)=15.00029
1297         zlev(3)=25.00176
1298         zlev(4)=35.00541
1299         zlev(5)=45.01332
1300         zlev(6)=55.0295
1301         zlev(7)=65.06181
1302         zlev(8)=75.12551
1303         zlev(9)=85.25037
1304         zlev(10)=95.49429
1305         zlev(11)=105.9699
1306         zlev(12)=116.8962
1307         zlev(13)=128.6979
1308         zlev(14)=142.1953
1309         zlev(15)=158.9606
1310         zlev(16)=181.9628
1311         zlev(17)=216.6479
1312         zlev(18)=272.4767
1313         zlev(19)=364.303
1314         zlev(20)=511.5348
1315         zlev(21)=732.2009
1316         zlev(22)=1033.217
1317         zlev(23)=1405.698
1318         zlev(24)=1830.885
1319         zlev(25)=2289.768
1320         zlev(26)=2768.242
1321         zlev(27)=3257.479
1322         zlev(28)=3752.442
1323         zlev(29)=4250.401
1324         zlev(30)=4749.913
1325         zlev(31)=5250.227
1326      ELSEIF (nmlev==1) THEN
1327         zlev(1)=0.0
1328      ELSE
1329         WRITE(*,*) 'Unknown number of levels'
1330         CALL abort
1331      ENDIF
1332
1333   END SUBROUTINE getlevs
1334
1335END PROGRAM fbstat
Note: See TracBrowser for help on using the repository browser.