Subversion Repositories lagranto.ecmwf

Rev

Rev 19 | Rev 23 | Go to most recent revision | Details | Compare with Previous | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
c     ************************************************************
2
c     * This package provides input routines to read the wind    *
3
c     * and other fields from IVE necdf files. The routines are  *
4
c     *                                                          *
5
c     * 1) input_open  : to open a data file                     *
6
c     * 2) input_grid  : to read the grid information, including *
7
c     *                  the vertical levels                     *
8
c     * 3) input_wind  : to read the wind components             *
9
c     * 4) input_close : to close an input file                  *
10
c     *                                                          *
11
c     * The file is characterised by an filename <filename> and  *
12
c     * a file identifier <fid>. The horizontal grid is given by *
13
c     * <xmin,xmax,ymin,ymax,dx,dy,nx,ny> where the pole of the  *
14
c     * rotated grid is given by <pollon,pollat>. The vertical   *
15
c     * grid is characterised by the surface pressure <ps> and   *
16
c     * the pressure at staggered <slev> and unstaggered <ulev>  *
17
c     * levels. The number of levels is given by <nz>. Finally,  *
18
c     * the retrieval of the wind <field> with name <fieldname>  *
19
c     * is characterised by a <time> and a missing data value    *
20
c     * <mdv>.                                                   *
21
c     *                                                          *
22
c     * Author: Michael Sprenger, Autumn 2008                    *
23
c     ************************************************************
24
 
25
c     ------------------------------------------------------------
26
c     Open input file
27
c     ------------------------------------------------------------
28
 
29
      subroutine input_open (fid,filename)
30
 
31
c     Open the input file with filename <filename> and return the
32
c     file identifier <fid> for further reference. 
33
 
21 michaesp 34
      use netcdf
35
 
3 michaesp 36
      implicit none
37
 
38
c     Declaration of subroutine parameters
39
      integer      fid              ! File identifier
40
      character*80 filename         ! Filename
41
 
42
c     Declaration of auxiliary variables
43
      integer      ierr
44
 
21 michaesp 45
c     Open netcdf file
46
      ierr = NF90_OPEN(TRIM(filename),nf90_nowrite, fid)
47
      IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
48
 
3 michaesp 49
c     Exception handling
50
      return
51
 
52
      end
21 michaesp 53
 
3 michaesp 54
c     ------------------------------------------------------------
55
c     Read information about the grid
56
c     ------------------------------------------------------------
57
 
58
      subroutine input_grid 
59
     >                   (fid,fieldname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
60
     >                    time,pollon,pollat,p3,ps,nz,ak,bk,stagz,
61
     >                    timecheck)
62
 
63
c     Read grid information at <time> from file with identifier <fid>. 
64
c     The horizontal grid is characterized by <xmin,xmax,ymin,ymax,dx,dy>
65
c     with pole position at <pollon,pollat> and grid dimension <nx,ny>.
66
c     The 3d arrays <p3(nx,ny,nz)> gives the vertical coordinates, either
67
c     on the staggered or unstaggered grid (with <stagz> as the flag).
68
c     The surface pressure is given in <ps(nx,ny)>. If <fid> is negative, 
69
c     only the grid dimensions and grid parameters (xmin...pollat,nz) are 
70
c     determined and returned (this is needed for dynamical allocation of 
71
c     memory).
72
 
21 michaesp 73
      use netcdf
74
 
3 michaesp 75
      implicit none
76
 
77
c     Declaration of subroutine parameters 
78
      integer      fid                 ! File identifier
79
      real         xmin,xmax,ymin,ymax ! Domain size
80
      real         dx,dy               ! Horizontal resolution
81
      integer      nx,ny,nz            ! Grid dimensions
82
      real         pollon,pollat       ! Longitude and latitude of pole
83
      real         p3(nx,ny,nz)        ! Staggered levels
84
      real         ps(nx,ny)           ! Surface pressure
85
      real         time                ! Time of the grid information
86
      real         ak(nz),bk(nz)       ! Ak and Bk for layers or levels
87
      real         stagz               ! Vertical staggering (0 or -0.5)
88
      character*80 fieldname           ! Variable from which to take grid info
89
      character*80 timecheck           ! Either 'yes' or 'no'
90
 
91
c     Numerical and physical parameters
92
      real          eps                 ! Numerical epsilon
93
      parameter    (eps=0.001)
94
 
95
c     Netcdf variables
96
      integer      vardim(4)
97
      real         varmin(4),varmax(4)
98
      real         mdv
99
      real         stag(4)
100
      integer      ndim
101
      character*80 cstfile
102
      integer      cstid
21 michaesp 103
      real         aklay(nz),bklay(nz)
3 michaesp 104
      integer      nvars
105
      character*80 vars(100)
21 michaesp 106
      integer        dimids (nf90_max_var_dims)
107
      character*80   dimname(nf90_max_var_dims)
108
      real,allocatable, dimension (:)     :: lon,lat,lev
109
      real,allocatable, dimension (:)     :: times
110
      real,allocatable, dimension (:,:)   :: tmp2
111
      real,allocatable, dimension (:,:,:) :: tmp3
112
      character*80  units
3 michaesp 113
 
21 michaesp 114
c     Auxiliary variables
3 michaesp 115
      integer      ierr       
116
      integer      i,j,k
117
      integer      isok
118
      real         tmp(200)
119
      character*80 varname
120
      real         rtime
21 michaesp 121
      integer      varid
122
      integer      cdfid
123
      integer      stat
124
      real         delta
125
      integer      closear
126
      real         maxps,minps
3 michaesp 127
 
21 michaesp 128
c     ------ Set file identifier --------------------------------------
3 michaesp 129
      if (fid.lt.0) then
21 michaesp 130
        cdfid = -fid
131
      else 
132
        cdfid = fid
133
      endif
3 michaesp 134
 
21 michaesp 135
c     Special handling if 3D pressure is
136
      if ( fieldname.eq.'P' ) then
137
         fieldname = 'U'
138
      endif
3 michaesp 139
 
21 michaesp 140
c     Get number of dimensions of variable -> ndim
141
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
142
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
143
      ierr = nf90_inquire_variable(cdfid, varid, ndims  = ndim)
144
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
145
      if ( ndim.ne.4 ) then
146
          print*,' ERROR: netCDF variables need to be 4D'
147
          print*,'      ',trim(fieldname)
148
          stop
149
      endif
3 michaesp 150
 
21 michaesp 151
c     Get dimensions -> vardim(1:ndim),dimname(1:ndim)
152
      ierr = NF90_INQ_VARID(cdfid,fieldname,varid)
153
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
154
      ierr = nf90_inquire_variable(cdfid, varid, 
155
     >                                   dimids = dimids(1:ndim))
156
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
157
      do i=1,ndim
158
           ierr = nf90_inquire_dimension(cdfid, dimids(i), 
159
     >                               name = dimname(i) )
160
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
161
           ierr = nf90_inquire_dimension(cdfid, dimids(i),len=vardim(i))
162
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
163
      enddo
3 michaesp 164
 
21 michaesp 165
c     Check whether the list of dimensions is OK
166
      if ( ( dimname(1).ne.'lon'  ).or.
167
     >     ( dimname(2).ne.'lat'  ).or. 
168
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
169
     >     ( dimname(4).ne.'time' ) )
170
     >then
171
        print*,' ERROR: the dimensions of the variable are not correct'
172
        print*,'        expected -> lon / lat / lev / time'
173
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
174
        stop
175
      endif
3 michaesp 176
 
21 michaesp 177
c     Allocate memory for reading arrays
178
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
179
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
180
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
181
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
182
      allocate(lon(vardim(1)),stat=stat)
183
      if (stat.ne.0) print*,'*** error allocating array lon     ***' 
184
      allocate(lat(vardim(2)),stat=stat)
185
      if (stat.ne.0) print*,'*** error allocating array lat     ***' 
186
      allocate(times(vardim(4)),stat=stat)
187
      if (stat.ne.0) print*,'*** error allocating array times   ***'
19 michaesp 188
 
21 michaesp 189
c     Get domain longitudes and latitudes
190
      varname = dimname(1)
191
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
192
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
193
      ierr = nf90_get_var(cdfid,varid,lon)
194
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
195
      varname = dimname(2)
196
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
197
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
198
      ierr = nf90_get_var(cdfid,varid,lat)
199
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
200
 
201
c     Get ak and bk
202
      varname='hyam'
203
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
204
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
205
      ierr = nf90_get_var(cdfid,varid,ak)
206
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
207
      varname='hybm'
208
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
209
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
210
      ierr = nf90_get_var(cdfid,varid,bk)
211
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
19 michaesp 212
 
21 michaesp 213
c     Check that unit of ak is in hPa - if necessary correct it
214
      varname='hyam'
215
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
216
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
217
      ierr = nf90_get_att(cdfid, varid, "units", units)
218
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
219
      if ( units.eq.'Pa' ) then
220
      print*,'Hallo'
221
         do i=1,vardim(3)
222
            ak(k) = 0.01 * ak(k)
223
         enddo
224
      endif
19 michaesp 225
 
21 michaesp 226
c     Get time information (check if time is correct)
227
      varname = 'time'
228
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
229
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
230
      ierr = nf90_get_var(cdfid,varid,times)
231
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
232
      isok=0
233
      do i=1,vardim(4)
234
        if (abs(time-times(i)).lt.eps) then
3 michaesp 235
               isok = 1
236
               rtime = times(i)
21 michaesp 237
        elseif (timecheck.eq.'no') then
3 michaesp 238
               isok = 1
239
               rtime = times(1)
21 michaesp 240
        endif
241
      enddo
242
      if ( isok.eq.0 ) then
243
         print*,' ERROR: time ',rtime,' not found on netCDF file' 
244
         stop
245
      endif
3 michaesp 246
 
21 michaesp 247
c     Read surface pressure
248
      varname='PS'
249
      ierr = NF90_INQ_VARID(cdfid,varname,varid)
250
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
251
      ierr = nf90_get_var(cdfid,varid,tmp2)
252
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
253
 
254
c     Check that surface pressure is in hPa
255
      maxps = -1.e39
256
      minps =  1.e39
257
      do i=1,vardim(1)
258
        do j=1,vardim(2)
259
             if (tmp2(i,j).gt.maxps) maxps = tmp2(i,j)
260
             if (tmp2(i,j).lt.minps) minps = tmp2(i,j)
261
        enddo
262
      enddo
263
      if ( (maxps.gt.1500.).or.(minps.lt.300.) ) then
264
         print*,' ERROR: surface pressre PS must be in hPa'
265
         print*,'       ',maxps,minps
266
         stop
267
      endif
268
 
269
c     Calculate layer and level pressures
270
      do i=1,vardim(1)
271
         do j=1,vardim(2)
272
               do k=1,vardim(3)
273
                  tmp3(i,j,k)=ak(k)+bk(k)*tmp2(i,j)
3 michaesp 274
               enddo
21 michaesp 275
         enddo
276
      enddo
3 michaesp 277
 
21 michaesp 278
c     Set the grid dimensions and constants
279
      nx      = vardim(1)
280
      ny      = vardim(2)
281
      nz      = vardim(3)
282
      xmin    = lon(1)
283
      ymin    = lat(1)
284
      xmax    = lon(nx)
285
      ymax    = lat(ny)
286
      dx      = (xmax-xmin)/real(nx-1)
287
      dy      = (ymax-ymin)/real(ny-1)
288
      pollon  = 0.
289
      pollat  = 90.
290
      stagz   = 0.
291
      delta   = xmax-xmin-360.
292
      if (abs(delta+dx).lt.eps) then
293
          xmax    = xmax + dx
294
          nx      = nx + 1
295
          closear = 1
296
      else
297
          closear = 0
298
      endif
3 michaesp 299
 
21 michaesp 300
c     Save the output arrays (if fid>0) - close arrays on request
301
      if ( fid.gt.0 ) then
3 michaesp 302
 
21 michaesp 303
         do j=1,vardim(2)
304
           do i=1,vardim(1)
305
             ps(i,j) = tmp2(i,j)
306
           enddo
307
           if (closear.eq.1) ps(vardim(1)+1,j) = ps(1,j)
3 michaesp 308
         enddo
309
 
21 michaesp 310
         do j=1,vardim(2)
311
           do k=1,vardim(3)
312
             do i=1,vardim(1)
313
               p3(i,j,k) = tmp3(i,j,k)
314
             enddo
315
             if (closear.eq.1) p3(vardim(1)+1,j,k) = p3(1,j,k)
316
           enddo
3 michaesp 317
         enddo
21 michaesp 318
      endif
3 michaesp 319
 
320
      return
321
 
322
      end
323
 
324
c     ------------------------------------------------------------
325
c     Read wind information
326
c     ------------------------------------------------------------
327
 
328
      subroutine input_wind (fid,fieldname,field,time,stagz,mdv,
329
     >                       xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,
330
     >                       timecheck)
331
 
332
c     Read the wind component <fieldname> from the file with identifier
333
c     <fid> and save it in the 3d array <field>. The vertical staggering 
334
c     information is provided in <stagz> and gives the reference to either
335
c     the layer or level field from <input_grid>. A consistency check is
336
c     performed to have an agreement with the grid specified by <xmin,xmax,
337
c     ymin,ymax,dx,dy,nx,ny,nz>.
338
 
21 michaesp 339
      use netcdf
340
 
3 michaesp 341
      implicit none
342
 
343
c     Declaration of variables and parameters
344
      integer      fid                 ! File identifier
345
      character*80 fieldname           ! Name of the wind field
346
      integer      nx,ny,nz            ! Dimension of fields
347
      real         field(nx,ny,nz)     ! 3d wind field
348
      real         stagz               ! Staggering in the z direction
349
      real         mdv                 ! Missing data flag
350
      real         xmin,xmax,ymin,ymax ! Domain size
351
      real         dx,dy               ! Horizontal resolution
352
      real         time                ! Time
353
      character*80 timecheck           ! Either 'yes' or 'no'
354
 
355
c     Numerical and physical parameters
356
      real        eps                 ! Numerical epsilon
357
      parameter  (eps=0.001)
358
      real        notimecheck         ! 'Flag' for no time check
359
      parameter  (notimecheck=7.26537)
360
 
361
c     Netcdf variables
362
      integer      ierr
363
      character*80 varname
364
      integer      vardim(4)
365
      real         varmin(4),varmax(4)
366
      real         stag(4)
367
      integer      ndim
368
      real         times(10)
369
      integer      ntimes
370
      character*80 cstfile
371
      integer      cstid
372
      real         aklay(200),bklay(200),aklev(200),bklev(200)
373
      real         ps(nx,ny)
21 michaesp 374
      integer      dimids (nf90_max_var_dims)
375
      character*80 dimname(nf90_max_var_dims)
376
      integer      varid
377
      integer      cdfid
378
      real,allocatable, dimension (:)     :: lon,lat
379
      real,allocatable, dimension (:,:)   :: tmp2
380
      real,allocatable, dimension (:,:,:) :: tmp3
3 michaesp 381
 
382
c     Auxiliary variables
383
      integer      isok
384
      integer      i,j,k
385
      integer      nz1
386
      real         rtime
21 michaesp 387
      integer      closear
388
      integer      stat
389
      real         delta
3 michaesp 390
 
21 michaesp 391
c     Get the number of dimensions -> ndim
392
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
393
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
394
      ierr = nf90_inquire_variable(fid, varid, ndims  = ndim)
395
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
19 michaesp 396
 
21 michaesp 397
c     Get the dimensions of the arrays -> varid(1:ndim)
398
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
399
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
400
      ierr = nf90_inquire_variable(fid, varid, 
401
     >                                   dimids = dimids(1:ndim))
402
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
403
      do i=1,ndim
404
           ierr = nf90_inquire_dimension(fid, dimids(i), 
405
     >                               name = dimname(i) )
406
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
407
           ierr = nf90_inquire_dimension(fid, dimids(i),len=vardim(i))
408
           IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
3 michaesp 409
      enddo
410
 
21 michaesp 411
c     Check whether the list of dimensions is OK
412
      if ( ( dimname(1).ne.'lon'  ).or.
413
     >     ( dimname(2).ne.'lat'  ).or.
414
     >     ( dimname(3).ne.'lev'  ).and.( dimname(3).ne.'lev_2'  ).or.
415
     >     ( dimname(4).ne.'time' ) )
416
     >then
417
        print*,' ERROR: the dimensions of the variable are not correct'
418
        print*,'        expected -> lon / lat / lev / time'
419
        print*, ( trim(dimname(i))//' / ',i=1,ndim )
420
        stop
421
      endif
3 michaesp 422
 
21 michaesp 423
c     Allocate memory for reading arrays - depending on <closear>
424
      allocate(tmp2(vardim(1),vardim(2)),stat=stat)
425
      if (stat.ne.0) print*,'*** error allocating array tmp2     ***'
426
      allocate(tmp3(vardim(1),vardim(2),vardim(3)),stat=stat)
427
      if (stat.ne.0) print*,'*** error allocating array tmp3     ***'
428
      allocate(lon(vardim(1)),stat=stat)
429
      if (stat.ne.0) print*,'*** error allocating array lon     ***'
430
      allocate(lat(vardim(2)),stat=stat)
431
      if (stat.ne.0) print*,'*** error allocating array lat     ***'
3 michaesp 432
 
21 michaesp 433
c     Get domain boundaries
434
      varname = dimname(1)
435
      ierr = NF90_INQ_VARID(fid,varname,varid)
436
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
437
      ierr = nf90_get_var(fid,varid,lon)
438
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
439
      varname = dimname(2)
440
      ierr = NF90_INQ_VARID(fid,varname,varid)
441
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
442
      ierr = nf90_get_var(fid,varid,lat)
443
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
3 michaesp 444
 
21 michaesp 445
c     Read data 
446
      ierr = NF90_INQ_VARID(fid,fieldname,varid)
447
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
448
      ierr = nf90_get_var(fid,varid,tmp3)
449
      IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
450
 
3 michaesp 451
c     If the field is 2D, expand it to 3D - simple handling of 2D tracing
452
      if ( vardim(3).eq.1 ) then
21 michaesp 453
         do i=1,vardim(1)
454
            do j=1,vardim(2)
455
               do k=1,vardim(3)
456
                  tmp3(i,j,k) = tmp3(i,j,1)
3 michaesp 457
               enddo
458
            enddo
459
         enddo
460
      endif
461
 
21 michaesp 462
c     Save the ouput array - close on request
463
      dx    = (varmax(1)-varmin(1))/real(vardim(1)-1)
464
      delta = varmax(1)-varmin(1)-360.
465
      if (abs(delta+dx).lt.eps) then
466
          closear = 1
467
      else
468
          closear = 0
469
      endif
3 michaesp 470
 
21 michaesp 471
      do j=1,vardim(2)
472
        do k=1,vardim(3)
473
          do i=1,vardim(1)
474
             field(i,j,k) = tmp3(i,j,k)
475
          enddo
476
          if (closear.eq.1) field(vardim(1)+1,j,k) = field(1,j,k)
477
        enddo
478
      enddo
479
 
480
c     Exit point
3 michaesp 481
      return
21 michaesp 482
 
3 michaesp 483
      end
484
 
485
c     ------------------------------------------------------------
486
c     Close input file
487
c     ------------------------------------------------------------
488
 
489
      subroutine input_close(fid)
490
 
491
c     Close the input file with file identifier <fid>.
492
 
21 michaesp 493
      use netcdf
494
 
3 michaesp 495
      implicit none
496
 
497
c     Declaration of subroutine parameters
498
      integer fid
499
 
500
c     Auxiliary variables
501
      integer ierr
502
 
503
c     Close file
21 michaesp 504
      ierr = NF90_CLOSE(fid)
505
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
506
 
507
      end
508
 
509
c     ------------------------------------------------------------
510
c     Get a list of variables on netCDF file
511
c     ------------------------------------------------------------
512
 
513
      subroutine input_getvars(fid,vnam,nvars)
514
 
515
c     List of variables on netCDF file
516
 
517
      use netcdf
518
 
519
      implicit none
520
 
521
c     Declaration of subroutine parameters
522
      integer      fid
523
      integer      nvars
524
      character*80 vnam(200)
525
 
526
c     Auxiliary variables
527
      integer ierr
528
      integer i
529
      integer nDims,nGlobalAtts,unlimdimid
530
 
531
      ierr = nf90_inquire(fid, nDims, nVars, nGlobalAtts, unlimdimid)
532
      IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
533
 
534
      do i=1,nVars
535
         ierr = nf90_inquire_variable(fid, i, name = vnam(i))
536
      enddo
3 michaesp 537
 
538
      end