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 branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/TOOLS/OBSTOOLS/src – NEMO

source: branches/2014/dev_r4650_UKMO12_CFL_diags_take2/NEMOGCM/TOOLS/OBSTOOLS/src/fbstat.F90 @ 5947

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

Reinstate svn Id keywords before merge

  • Property svn:executable set to *
  • Property svn:keywords set to Id
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.