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/dev_2802_OBStools/NEMOGCM/TOOLS/OBSTOOLS – NEMO

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

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

Adding obs tools to branch
  • Property svn:executable set to *
File size: 20.5 KB
Line 
1PROGRAM fbstat
2   USE obs_fbm
3   USE fbaccdata
4   USE coords
5   USE omonainfo
6   IMPLICIT NONE
7   TYPE(obfbdata) :: fbdata
8   CHARACTER(len=256) :: filename
9   INTEGER :: jfile,jbox,jlev,jfirst,jvar,jadd,ji
10#ifndef NOIARGCPROTO
11   INTEGER,EXTERNAL :: iargc
12#endif
13   REAL,DIMENSION(:),ALLOCATABLE :: zlev
14   INTEGER :: nmlev, nfiles
15   LOGICAL :: lexists,lomona,ltext
16   REAL, ALLOCATABLE, DIMENSION(:,:,:) :: zdat3d
17   REAL, ALLOCATABLE, DIMENSION(:,:) :: zdat2d
18   INTEGER,DIMENSION(1) :: itime
19   INTEGER :: inidate,icurdate,loopno,ityp,iloopno
20   INTEGER :: nvar,nadd
21   CHARACTER(len=4) :: expver
22   CHARACTER(len=7) :: cltyp
23   CHARACTER(len=128) :: cdfmthead,cdfmtbody
24   LOGICAL :: lnear,linner,linnerp,linnerini,lpassive,lhistogram,lfound
25   INTEGER :: nqc
26   CHARACTER(len=ilenname), DIMENSION(:), ALLOCATABLE :: cname,caddname
27   CHARACTER(len=20) :: carea
28   INTEGER :: jstabox, jendbox
29   NAMELIST/namfbstat/ltext,lomona,nmlev,inidate,icurdate,loopno,&
30      &               expver,lnear,linner,lpassive,lhistogram, &
31      &               zhistmax,zhistmin,zhiststep,zcheck,carea
32
33   ltext=.TRUE.
34   lomona=.FALSE.
35   nmlev=31
36   inidate=19010101
37   icurdate=19010116
38   loopno=-1
39   expver='test'
40   lnear=.TRUE.
41   linner=.FALSE.
42   lpassive=.FALSE.
43   lhistogram=.FALSE.
44   zhistmin(:)=-10.0
45   zhistmax(:)=10.0
46   zhiststep(:)=0.1
47   zcheck(:)=10.0
48   nqc=2
49   carea='all'
50   INQUIRE(file='namfbstat.in',exist=lexists)
51   IF (lexists) THEN
52      OPEN(10,file='namfbstat.in')
53      READ(10,namfbstat)
54      CLOSE(10)
55      WRITE(*,namfbstat)
56   ENDIF
57   IF (iargc()==0) THEN
58      WRITE(*,*)'Usage:'
59      WRITE(*,*)'fbstat [-nmlev] <filenames>'
60      CALL abort
61   ENDIF
62   jfirst=1
63   DO ji=1,2
64      CALL getarg(jfirst,filename)
65      IF (filename=='-42') THEN
66         nmlev=42
67         jfirst=jfirst+1
68      ELSEIF(filename=='-31') THEN
69         nmlev=31
70         jfirst=jfirst+1
71      ELSEIF(filename=='-1') THEN
72         nmlev=1
73         lnear=.TRUE.
74         jfirst=jfirst+1
75      ELSEIF(filename=='-q') THEN
76         jfirst=jfirst+1
77         CALL getarg(jfirst,filename)
78         READ(filename,'(I4)')nqc
79         IF ((nqc<0).OR.(nqc>4)) THEN
80            WRITE(*,*)'Quality control option (-q) should be 1 to 4'
81            CALL abort
82         ENDIF
83         jfirst=jfirst+1
84      ENDIF
85   END DO
86   nfiles=iargc()
87   
88   IF (carea/='all') THEN
89      IF (lomona) THEN
90         WRITE(*,*)'For omona files carea has to be all'
91         CALL abort
92      ENDIF
93      lfound=.FALSE.
94      DO jbox=1,nbox
95         IF (TRIM(carea)==TRIM(cl_boxes(jbox))) THEN
96            jstabox=jbox
97            jendbox=jbox
98            lfound=.TRUE.
99         ENDIF
100      ENDDO
101      IF (.NOT.lfound) THEN
102         WRITE(*,*)'Area not found'
103         CALL abort
104      ENDIF
105   ELSE
106      jstabox=1
107      jendbox=nbox
108   ENDIF
109   PRINT *,'jstabox,jendbox=',jstabox,jendbox
110
111   IF (.NOT.lnear) nmlev=nmlev-1
112
113   ALLOCATE(&
114      & zlev(nmlev) &
115      & )
116   IF(lnear) THEN
117      CALL getlevs(nmlev,zlev)
118   ELSE
119      CALL getlevsmean(nmlev,zlev)
120   ENDIF
121
122   DO jfile=jfirst, nfiles
123      CALL getarg(jfile,filename)
124      WRITE(*,*)'Handling file : ',TRIM(filename)
125      CALL read_obfbdata(TRIM(filename),fbdata)
126      IF (jfile==jfirst) THEN
127         nvar=fbdata%nvar
128         nadd=fbdata%nadd
129         IF (lhistogram) THEN
130            IF (nvar>maxvars) THEN
131               WRITE(*,*)'fbstat.F90: Increase maxvars to ',nvar
132               WRITE(*,*)'if you want histograms'
133               CALL abort
134            ENDIF
135            DO jvar = 1, nvar
136               hist(jvar)%npoints=(zhistmax(jvar)-zhistmin(jvar))&
137                  &              /zhiststep(jvar)+1
138               WRITE(*,*)'Number of points in histogram      = ',&
139                  &  hist(jvar)%npoints
140               WRITE(*,*)'Size of histogram array (elements) = ',&
141                  &  hist(jvar)%npoints*nmlev*nbox*nadd
142               ALLOCATE(&
143                  & hist(jvar)%nhist(hist(jvar)%npoints,nmlev,nbox,nadd) &
144                  & )
145               hist(jvar)%nhist(:,:,:,:)=0
146            ENDDO
147         ENDIF
148         ALLOCATE(&
149            & inum(nmlev,nbox,nadd,nvar),  &
150            & zmean(nmlev,nbox,nadd,nvar), &
151            & zrms(nmlev,nbox,nadd,nvar),  &
152            & zsdev(nmlev,nbox,nadd,nvar), &
153            & cname(nvar),                 &
154            & caddname(nadd)               &
155            & )
156         DO jvar = 1, nvar
157            cname(jvar) = fbdata%cname(jvar)
158         END DO
159         DO jadd = 1, nadd
160            caddname(jadd) = fbdata%caddname(jadd)
161         END DO
162         inum(:,:,:,:)=0
163         zmean(:,:,:,:)=0.0
164         zrms(:,:,:,:)=0.0
165         zsdev(:,:,:,:)=0.0
166      ELSE
167         IF (fbdata%nvar/=nvar) THEN
168            WRITE(*,*)'Different number of nvar ',fbdata%nvar,' in ',&
169               & TRIM(filename)
170            CALL abort
171         ENDIF
172         IF (fbdata%nadd/=nadd) THEN
173            WRITE(*,*)'Different number of nadd ',fbdata%nadd,' in ',&
174               & TRIM(filename)
175            CALL abort
176         ENDIF
177      ENDIF
178      CALL fb_stat(fbdata,jstabox,jendbox,nmlev,zlev,lnear,nqc,lhistogram)
179      CALL dealloc_obfbdata(fbdata)
180   ENDDO
181
182   DO jvar=1, nvar
183      DO jadd=1, nadd
184         DO jbox=jstabox, jendbox
185            DO jlev = 1, nmlev
186               IF ( inum(jlev,jbox,jadd,jvar) > 0 ) THEN
187                  zmean(jlev,jbox,jadd,jvar) = &
188                     & zmean(jlev,jbox,jadd,jvar)/inum(jlev,jbox,jadd,jvar)
189                  zrms(jlev,jbox,jadd,jvar) = &
190                     & SQRT(zrms(jlev,jbox,jadd,jvar)/inum(jlev,jbox,jadd,jvar))
191                  zsdev(jlev,jbox,jadd,jvar) = &
192                     & SQRT(MAX(zrms(jlev,jbox,jadd,jvar)**2-zmean(jlev,jbox,jadd,jvar)**2,0.0))
193               ELSE
194                  zmean(jlev,jbox,jadd,jvar) = fbrmdi
195                  zrms(jlev,jbox,jadd,jvar) = fbrmdi
196                  zsdev(jlev,jbox,jadd,jvar) = fbrmdi
197               ENDIF
198            ENDDO
199         ENDDO
200      ENDDO
201   ENDDO
202
203   IF (ltext) THEN
204      DO jvar=1, nvar
205         DO jadd=1, nadd
206            DO jbox=jstabox, jendbox
207               WRITE(filename,'(5A)')TRIM(cname(jvar)),&
208                  &                  TRIM(caddname(jadd)),'_',&
209                  &                  TRIM(cl_boxes(jbox)),'.dat'
210               OPEN(10,file=TRIM(filename))
211               DO jlev = 1, nmlev
212                  WRITE(10,'(F16.7,2I12,3F17.10)') zlev(jlev), &
213                     & jlev, inum(jlev,jbox,jadd,jvar), &
214                     & zmean(jlev,jbox,jadd,jvar), &
215                     & zrms(jlev,jbox,jadd,jvar), &
216                     & zsdev(jlev,jbox,jadd,jvar)
217               ENDDO
218               CLOSE(10)
219            ENDDO
220         ENDDO
221      ENDDO
222   ENDIF
223
224   IF (lhistogram) THEN
225      DO jvar=1, nvar
226         DO jadd=1, nadd
227            DO jbox=jstabox, jendbox
228               WRITE(filename,'(5A)')TRIM(cname(jvar)),&
229                  &                  TRIM(caddname(jadd)),'_',&
230                  &                  TRIM(cl_boxes(jbox)),'_histogram.dat'
231               OPEN(10,file=TRIM(filename))
232               WRITE(10,'(A10,100F10.2)')'#',(zlev(jlev),jlev=1,nmlev)
233               DO ji=1,hist(jvar)%npoints
234                  WRITE(10,'(F10.2,100I10)') &
235                     & zhistmin(jvar)+(ji-1)*zhiststep(jvar), &
236                     & (hist(jvar)%nhist(ji,jlev,jbox,jadd),jlev=1,nmlev)
237               ENDDO
238               CLOSE(10)
239            ENDDO
240         ENDDO
241      ENDDO
242   ENDIF
243
244   IF (lomona) THEN
245
246      IF (nmlev>1) THEN
247         ALLOCATE(zdat3d(nmlev,nbox,1))
248      ELSE
249         ALLOCATE(zdat2d(nbox,1))
250      ENDIF
251
252      cl_expnam=expver
253      WRITE(cl_date,'(I8.8)')inidate
254      i_dp = nmlev
255      itime=icurdate
256      linnerp=.TRUE.
257      iloopno = loopno
258      linnerini = linner
259      DO jvar = 1, nvar
260         linner = linnerini
261         loopno = iloopno
262         SELECT CASE (TRIM(cname(jvar)))
263         CASE('POTM')
264            cl_var = 'votemper'
265         CASE('PSAL')
266            cl_var='vosaline'
267         CASE('SLA')
268            cl_var='sossheig'
269         CASE('SST')
270            cl_var='sosstsst'
271         END SELECT
272         DO jadd = 1, nadd
273            linner = (caddname(jadd)(1:3)=='Hxa').OR.linner
274            IF (lpassive) THEN
275               ityp=145
276            ELSE
277               IF (linner) THEN
278                  linnerp=.TRUE.
279                  ityp=144
280                  IF (jadd>1) loopno=loopno+1
281               ELSE
282                  ityp=142
283                  IF (.NOT.linnerp) THEN
284                     IF (jadd>1) loopno=loopno+1
285                  ENDIF
286               ENDIF
287            ENDIF
288            WRITE(cltyp,'(I3.3,A1,I2.2)')ityp,'_',loopno
289            CALL obs_variable_att(cltyp)
290            IF (nmlev>1) THEN
291               zdat3d(:,:,1) = zmean(:,:,jadd,jvar)
292               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
293                  &                    cl_boxes,REAL(fbrmdi))
294               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
295            ELSE
296               zdat2d(:,1) = zmean(1,:,jadd,jvar)
297               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
298                  &                    cl_boxes,REAL(fbrmdi))
299            ENDIF
300            IF (lpassive) THEN
301               ityp=245
302            ELSE
303               IF (linner) THEN
304                  ityp=244
305               ELSE
306                  ityp=242
307               ENDIF
308            ENDIF
309            WRITE(cltyp,'(I3.3,A1,I2.2)')ityp,'_',loopno
310            CALL obs_variable_att(cltyp)
311            IF (nmlev>1) THEN
312               zdat3d(:,:,1) = zrms(:,:,jadd,jvar)
313               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
314                  &                    cl_boxes,REAL(fbrmdi))
315               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
316            ELSE
317               zdat2d(:,1) = zrms(1,:,jadd,jvar)
318               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
319                  &                    cl_boxes,REAL(fbrmdi))
320            ENDIF
321            IF (lpassive) THEN
322               ityp=345
323            ELSE
324               IF (linner) THEN
325                  ityp=344
326               ELSE
327                  ityp=342
328               ENDIF
329            ENDIF
330            WRITE(cltyp,'(I3.3,A1,I2.2)')ityp,'_',loopno
331            CALL obs_variable_att(cltyp)
332            IF (nmlev>1) THEN
333               zdat3d(:,:,1) = zsdev(:,:,jadd,jvar)
334               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
335                  &                    cl_boxes,REAL(fbrmdi))
336               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
337            ELSE
338               zdat2d(:,1) = zsdev(1,:,jadd,jvar)
339               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
340                  &                    cl_boxes,REAL(fbrmdi))
341            ENDIF
342            IF (lpassive) THEN
343               ityp=445
344            ELSE
345               IF (linner) THEN
346                  ityp=444
347               ELSE
348                  ityp=442
349               ENDIF
350            ENDIF
351            WRITE(cltyp,'(I3.3,A1,I2.2)')ityp,'_',loopno
352            CALL obs_variable_att(cltyp)
353            IF (nmlev>1) THEN
354               zdat3d(:,:,1) = inum(:,:,jadd,jvar)
355               CALL write_omona_netcdf(cl_filename_out,zdat3d,itime, &
356                  &                    cl_boxes,REAL(fbrmdi))
357               CALL write_dep_netcdf(cl_filename_out,cl_boxes,zlev)
358            ELSE
359               zdat2d(:,1) = inum(1,:,jadd,jvar)
360               CALL write_omona_netcdf(cl_filename_out,zdat2d,itime, &
361                  &                    cl_boxes,REAL(fbrmdi))
362            ENDIF
363            linner=.FALSE.
364         ENDDO
365      ENDDO
366
367      IF (nmlev>1) THEN
368         DEALLOCATE(zdat3d)
369      ELSE
370         DEALLOCATE(zdat2d)
371      ENDIF
372     
373   ENDIF
374
375CONTAINS
376
377   SUBROUTINE fb_stat(fbdata,nstabox, nendbox, nmlev,zlev,lnear,kqc,lhist)
378      USE fbaccdata
379      USE coords
380      TYPE(obfbdata) :: fbdata
381      INTEGER :: nstabox,nendbox
382      INTEGER :: nmlev
383      REAL :: zlev(nmlev)
384      LOGICAL :: lnear
385      INTEGER :: kqc
386      LOGICAL :: lhist
387      INTEGER :: jlev, jobs, jvar, klev,jlev2,ih
388      REAL :: zarea(4),zlat,zlon,zdiff,zdiff2
389     
390      DO jbox = nstabox, nendbox
391         CALL coord_area(cl_boxes(jbox),zarea)
392         DO jobs = 1, fbdata%nobs
393            zlat = fbdata%pphi(jobs)
394            zlon = fbdata%plam(jobs)
395            IF (zlon<0) zlon=zlon+360
396            IF (zlon>360) zlon=zlon-360
397            IF ( ( zlat .GE. zarea(3) ) .AND. &
398               & ( zlat .LE. zarea(4) ) .AND. &
399               & ( ( ( zlon .GE. zarea(1) ) .AND. &
400               &     ( zlon .LE. zarea(2) ) ) .OR. &
401               &   ( ( zarea(2) .LE. zarea(1) ) .AND. &
402               &     ( zlon .GE. zarea(1) ) .AND. &
403               &     ( zlon .LE. 360        ) ) .OR. &
404               &   ( ( zarea(2)   .LE. zarea(1) ) .AND. &
405               &     ( zlon .GE. 0          ) .AND. &
406               &     ( zlon .LE. zarea(2) ) ) ) ) THEN
407           
408               DO jlev = 1, fbdata%nlev
409                  DO jvar = 1, fbdata%nvar
410                     IF (nmlev==1) THEN
411                        klev=1
412                     ELSE
413                        IF (lnear) THEN
414                           zdiff=ABS(fbdata%pdep(jlev,jobs)-zlev(1))
415                           klev=1
416                           DO jlev2=2,nmlev
417                              zdiff2=ABS(fbdata%pdep(jlev,jobs)-zlev(jlev2))
418                              IF (zdiff2<zdiff) THEN
419                                 klev=jlev2
420                                 zdiff=zdiff2
421                              ENDIF
422                           ENDDO
423                        ELSE
424                           klev = fbdata%iobsk(jlev,jobs,jvar)-1
425                        ENDIF
426                        IF ( klev > nmlev ) THEN
427                           DO jadd = 1, fbdata%nvar
428                              IF ( ABS(fbdata%padd(jlev,jobs,jadd,jvar))<9000 ) THEN
429                                 WRITE(*,*)'Error in fb_stat'
430                                 WRITE(*,*)'Increase nmlev to at least ',klev
431                                 klev=nmlev
432                                 CALL abort
433                              ENDIF
434                           ENDDO
435                        ENDIF
436                     ENDIF
437                     IF ( fbdata%ivlqc(jlev,jobs,jvar) > kqc ) CYCLE
438                     IF (( klev > 0 ).AND. &
439                        &(ABS(fbdata%pob(jlev,jobs,jvar)) < 9000 )) THEN
440                        DO jadd = 1, fbdata%nadd
441                           IF ( ABS(fbdata%padd(jlev,jobs,jadd,jvar)) < 9000 ) THEN
442                              zdiff = ( fbdata%padd(jlev,jobs,jadd,jvar) - &
443                                 &      fbdata%pob(jlev,jobs,jvar) )
444                              inum(klev,jbox,jadd,jvar) = inum(klev,jbox,jadd,jvar) + 1
445                              zmean(klev,jbox,jadd,jvar) = zmean(klev,jbox,jadd,jvar) + &
446                                 & zdiff
447                              zrms(klev,jbox,jadd,jvar) = zrms(klev,jbox,jadd,jvar) + &
448                                 & zdiff * zdiff
449                              IF (ABS(zdiff)>zcheck(jvar)) THEN
450                                 WRITE(*,*)'Departure outside check range ',&
451                                    & TRIM(fbdata%cname(jvar)),' entry ',&
452                                    & fbdata%caddname(jadd)
453                                 WRITE(*,*)'Depar = ',zdiff
454                                 WRITE(*,*)'Check = ',zcheck(jvar)
455                                 WRITE(*,*)'Id    = ',fbdata%cdwmo(jobs)
456                                 WRITE(*,*)'Lat   = ',fbdata%pphi(jobs)
457                                 WRITE(*,*)'Lon   = ',fbdata%plam(jobs)
458                                 WRITE(*,*)'Tim   = ',fbdata%ptim(jobs)
459                                 WRITE(*,*)'Depth = ',fbdata%pdep(jlev,jobs)
460                                 WRITE(*,*)'Obs   = ',fbdata%pob(jlev,jobs,jvar)
461                                 WRITE(*,*)'Var   = ',fbdata%padd(jlev,jobs,jadd,jvar)
462                                 WRITE(*,*)'QC    = ',fbdata%ivlqc(jlev,jobs,jvar)
463                                 WRITE(*,*)'QCflag= ',fbdata%ivlqcf(:,jlev,jobs,jvar)
464                              ENDIF
465                              IF (lhist) THEN
466                                 ih=NINT((zdiff-zhistmin(jvar))/zhiststep(jvar))+1
467                                 IF ((ih>=1).AND.(ih<=hist(jvar)%npoints)) THEN
468                                    hist(jvar)%nhist(ih,klev,jbox,jadd) = &
469                                       hist(jvar)%nhist(ih,klev,jbox,jadd) +1
470                                 ELSE
471                                    WRITE(*,*)'Histogram value outside range for ',&
472                                       & TRIM(fbdata%cname(jvar)),' entry ',&
473                                       & fbdata%caddname(jadd)
474                                    WRITE(*,*)'Value = ',zdiff
475                                    WRITE(*,*)'Range = ',zhistmin(jvar),zhistmax(jvar)
476                                    WRITE(*,*)'Step  = ',zhiststep(jvar)
477                                    WRITE(*,*)'Index = ',ih
478                                    WRITE(*,*)'Id    = ',fbdata%cdwmo(jobs)
479                                    WRITE(*,*)'Lat   = ',fbdata%pphi(jobs)
480                                    WRITE(*,*)'Lon   = ',fbdata%plam(jobs)
481                                    WRITE(*,*)'Tim   = ',fbdata%ptim(jobs)
482                                    WRITE(*,*)'Depth = ',fbdata%pdep(jlev,jobs)
483                                    WRITE(*,*)'Obs   = ',fbdata%pob(jlev,jobs,jvar)
484                                    WRITE(*,*)'Var   = ',fbdata%padd(jlev,jobs,jadd,jvar)
485                                    WRITE(*,*)'QC    = ',fbdata%ivlqc(jlev,jobs,jvar)
486                                    WRITE(*,*)'QCflag= ',fbdata%ivlqcf(:,jlev,jobs,jvar)
487                                 ENDIF
488                              ENDIF
489                           ENDIF
490                        ENDDO
491                     ENDIF
492                  ENDDO
493               ENDDO
494            ENDIF
495         ENDDO
496      ENDDO
497     
498   END SUBROUTINE fb_stat
499
500   SUBROUTINE getlevsmean(nmlev,zlev)
501      IMPLICIT NONE
502      INTEGER :: nmlev
503      REAL,DIMENSION(nmlev) :: zlev
504      REAL,DIMENSION(nmlev+1) :: ztmp
505      INTEGER :: i
506
507      zlev(:)=9999.9
508      CALL getlevs(nmlev+1,ztmp)
509      DO i=1,nmlev
510         zlev(i)=0.5*(ztmp(i)+ztmp(i+1))
511      ENDDO
512
513   END SUBROUTINE getlevsmean
514
515   SUBROUTINE getlevs(nmlev,zlev)
516      IMPLICIT NONE
517      INTEGER :: nmlev
518      REAL,DIMENSION(nmlev) :: zlev
519     
520      zlev(:)=9999.9
521
522      IF (nmlev==42) THEN
523         zlev(1)=5.02159
524         zlev(2)=15.07854
525         zlev(3)=25.16046
526         zlev(4)=35.27829
527         zlev(5)=45.44776
528         zlev(6)=55.69149
529         zlev(7)=66.04198
530         zlev(8)=76.54591
531         zlev(9)=87.27029
532         zlev(10)=98.31118
533         zlev(11)=109.8062
534         zlev(12)=121.9519
535         zlev(13)=135.0285
536         zlev(14)=149.4337
537         zlev(15)=165.7285
538         zlev(16)=184.6975
539         zlev(17)=207.4254
540         zlev(18)=235.3862
541         zlev(19)=270.5341
542         zlev(20)=315.3741
543         zlev(21)=372.9655
544         zlev(22)=446.8009
545         zlev(23)=540.5022
546         zlev(24)=657.3229
547         zlev(25)=799.5496
548         zlev(26)=967.9958
549         zlev(27)=1161.806
550         zlev(28)=1378.661
551         zlev(29)=1615.291
552         zlev(30)=1868.071
553         zlev(31)=2133.517
554         zlev(32)=2408.583
555         zlev(33)=2690.780
556         zlev(34)=2978.166
557         zlev(35)=3269.278
558         zlev(36)=3563.041
559         zlev(37)=3858.676
560         zlev(38)=4155.628
561         zlev(39)=4453.502
562         zlev(40)=4752.021
563         zlev(41)=5050.990
564         zlev(42)=5350.272
565      ELSEIF (nmlev==31) THEN
566         zlev(1)=4.999938
567         zlev(2)=15.00029
568         zlev(3)=25.00176
569         zlev(4)=35.00541
570         zlev(5)=45.01332
571         zlev(6)=55.0295
572         zlev(7)=65.06181
573         zlev(8)=75.12551
574         zlev(9)=85.25037
575         zlev(10)=95.49429
576         zlev(11)=105.9699
577         zlev(12)=116.8962
578         zlev(13)=128.6979
579         zlev(14)=142.1953
580         zlev(15)=158.9606
581         zlev(16)=181.9628
582         zlev(17)=216.6479
583         zlev(18)=272.4767
584         zlev(19)=364.303
585         zlev(20)=511.5348
586         zlev(21)=732.2009
587         zlev(22)=1033.217
588         zlev(23)=1405.698
589         zlev(24)=1830.885
590         zlev(25)=2289.768
591         zlev(26)=2768.242
592         zlev(27)=3257.479
593         zlev(28)=3752.442
594         zlev(29)=4250.401
595         zlev(30)=4749.913
596         zlev(31)=5250.227
597      ELSEIF (nmlev==1) THEN
598         zlev(1)=0.0
599      ELSE
600         WRITE(*,*) 'Unknown number of levels'
601         CALL abort
602      ENDIF
603
604   END SUBROUTINE getlevs
605
606END PROGRAM fbstat
Note: See TracBrowser for help on using the repository browser.