Subversion Repositories lagranto.20cr

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
4 michaesp 1
      PROGRAM density
2
 
3
      use netcdf
4
 
5
      implicit none
6
 
7
c     ---------------------------------------------------------------------
8
c     Declaration of variables
9
c     ---------------------------------------------------------------------
10
 
11
c     Parameter and working arrays
12
      real                                    radius
13
      character*80                            runit
14
      integer                                 nx,ny
15
      integer                                 nlonlat
16
      real                                    dlonlat
17
      real                                    xmin,ymin,dx,dy
18
      real                                    clon,clat
19
      integer                                 ntime,nfield,ntra
20
      character*80                            inpfile
21
      character*80                            outfile
22
      character*80                            mode
23
      real                                    param
24
      integer                                 opts,npts
25
      integer                                 step
26
      character*80                            gridtype
27
      character*80                            field
28
      integer                                 crefile,crevar
29
      real,allocatable,    dimension (:,:) :: cnt,res,fld,area
30
      real,allocatable,    dimension (:)   :: traj
31
      real,allocatable,    dimension (:)   :: olon,olat,otim,ofld
32
      real,allocatable,    dimension (:)   :: nlon,nlat,ntim,nfld
33
 
34
c     Output format
35
      character*80                            outformat
36
 
37
c     Physical and mathematical constants
38
      real                                    pi180
39
      parameter                               (pi180=3.14159/180.)
40
      real                                    deltay
41
      parameter                               (deltay=111.)
42
      real                                    eps
43
      parameter                               (eps=0.001)
44
 
45
c     Input trajectories (see iotra.f)
46
      integer                                 inpmode
47
      real,allocatable, dimension (:,:,:) ::  trainp     
48
      integer                                 reftime(6)      
49
      character*80                            varsinp(100)   
50
      integer,allocatable, dimension (:) ::   sel_flag
51
      character*80                            sel_file
52
      character*80                            sel_format
53
 
54
c     Auxiliary variables
55
      character*80                            cdfname,varname
56
      integer                                 i,j,k
57
      integer                                 stat
58
      integer,allocatable, dimension (:,:) :: connect0
59
      integer                                 connectval0
60
      integer,allocatable, dimension (:,:) :: connect1
61
      integer                                 connectval1
62
      integer,allocatable, dimension (:,:) :: connect2
63
      integer                                 connectval2
64
      real                                    slat
65
      integer                                 ipre
66
      real                                    addvalue
67
      real                                    xmax,ymax
68
      real ,allocatable, dimension (:)  ::    odist,ndist
69
      real                                    dt
70
      integer                                 fid
71
      integer                                 dynamic_grid
72
      real                                    ycen,xcen
73
      integer                                 indx,indy
74
      character*80                            unit
75
      real                                    pollon,pollat
76
      real                                    rlon0,rlat0,rlon,rlat
77
      real                                    lon,lat
78
      real                                    crot
79
      integer                                 count
80
      character*80                            longname, varunit
81
      real                                    time
82
      integer                                 ind
83
      integer                                 ifield
84
      real                                    hhmm,frac
85
      integer                                 ierr,ncID
86
 
87
c     External functions
88
      real         lmstolm,lmtolms
89
      real         phstoph,phtophs
90
      external     lmstolm,lmtolms,phstoph,phtophs
91
 
92
      real         sdis
93
      external     sdis
94
 
95
c     ---------------------------------------------------------------------
96
c     Preparations
97
c     ---------------------------------------------------------------------
98
 
99
c     Write start message
100
      print*,'========================================================='
101
      print*,'              *** START OF PROGRAM DENSITY ***'
102
      print*
103
 
104
c     Read input parameters
105
      open(10,file='density.param')
106
       read(10,*) inpfile
107
       read(10,*) outfile
108
       read(10,*) field
109
       read(10,*) ntime,nfield,ntra
110
       read(10,*) gridtype
111
       if ( gridtype.eq.'latlon' ) then 
112
          read(10,*) nx,ny,xmin,ymin,dx,dy
113
       elseif ( gridtype.eq.'rotated') then
114
          read(10,*) clon,clat,nlonlat,dlonlat
115
       else
116
          print*,' ERROR: unsupported grid type ',trim(gridtype)
117
          stop
118
       endif
119
       read(10,*) radius,runit
120
       read(10,*) mode
121
       read(10,*) param
122
       read(10,*) step
123
       read(10,*) sel_file
124
       read(10,*) sel_format
125
       read(10,*) crefile
126
       read(10,*) crevar
127
      close(10)
128
 
129
c     Get the grid parameters if <crefile=0>
130
      if ( crefile.eq.0 ) then
131
 
132
           ierr = nf90_open  (trim(outfile), NF90_NOWRITE  , ncID)
133
 
134
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'grid'   ,gridtype ) 
135
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clon'   ,clon     )
136
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'clat'   ,clat     )
137
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nlonlat',nlonlat  )
138
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dlonlat',dlonlat  )
139
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'nx'     ,nx       )
140
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ny'     ,ny       )
141
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dx'     ,dx       )
142
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'dy'     ,dy       )
143
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'xmin'   ,xmin     )
144
           ierr = nf90_get_att(ncID, NF90_GLOBAL, 'ymin'   ,ymin     )
145
 
146
           ierr = nf90_close(ncID)
147
 
148
           print*,'**** GRID PARAMETERS IMPORTED ',
149
     >            'FROM NETCDF FILE!!!! ****'
150
           print*
151
 
152
      endif
153
 
154
c     Check for consistency
155
      if ( (step.ne.0).and.(mode.ne.'keep') ) then
156
         print*," ERROR: interpolation is only possible for all",
157
     >                   ' time steps... Stop'
158
         stop
159
      endif
160
 
161
c     Set the number of times (just code aesthetics)
162
      opts=ntime
163
 
164
c     Set grid parameters for rotated grid
165
      if ( gridtype.eq.'rotated' ) then
166
         nx   = nlonlat
167
         ny   = nlonlat
168
         dx   = dlonlat
169
         dy   = dlonlat
170
         xmin = - real(nlonlat-1)/2. * dx
171
         xmax = + real(nlonlat-1)/2. * dx
172
         ymin = - real(nlonlat-1)/2. * dy
173
         ymax = + real(nlonlat-1)/2. * dy
174
      endif
175
 
176
c     Set the flag for dynamic grid adjustment
177
      if ( (nx.eq.0).or.(ny.eq.0) ) then
178
         dynamic_grid = 1
179
      else
180
         dynamic_grid = 0
181
      endif
182
 
183
c     Print status information
184
      print*,'---- INPUT PARAMETERS -----------------------------------'
185
      print* 
186
      print*,'Input                : ',trim(inpfile)
187
      print*,'Output               : ',trim(outfile)
188
      print*,'Field                : ',trim(field)
189
      print*,'Trajectory           : ',ntime,nfield,ntra
190
      print*,'Grid type            : ',trim(gridtype)
191
      if ( dynamic_grid.eq.1 ) then
192
         print*,'Grid                 : dynamic (see below)'
193
      elseif ( gridtype.eq.'latlon' ) then
194
         print*,'Grid   nlon,nlat     : ',nx,ny
195
         print*,'       lonmin,latmin : ',xmin,ymin
196
         print*,'       dlon,dlat     : ',dx,dy
197
      elseif ( gridtype.eq.'rotated' ) then
198
         print*,'Grid   clon,clat     : ',clon,clat
199
         print*,'       nlonlat       : ',nlonlat
200
         print*,'       dlonlat       : ',dlonlat
201
      endif
202
      print*,'Filter radius        : ',radius,' ',trim(runit)
203
      print*,'Mode                 : ',trim(mode)
204
      if ( ( mode.eq.'time'  ).or.
205
     >     ( mode.eq.'space' ).or.
206
     >     (mode.eq.'grid' ) ) 
207
     >then
208
         print*,'Parameter            : ',param
209
      endif
210
      if ( step.eq.0 ) then
211
         print*,'Time step            : all'
212
      elseif (step.gt.0) then
213
         print*,'Time step            : ',step
214
      endif
215
      print*,'Selection file       : ',trim(sel_file)
216
      print*,'Selection format     : ',trim(sel_file)
217
      print*,'Flag <crefile>       : ',crefile
218
      print*,'Flag <crevar>        : ',crevar
219
 
220
c     Check whether mode is valid
221
      if ((mode.ne.'keep'  ).and.
222
     >    (mode.ne.'time'  ).and.
223
     >    (mode.ne.'space' ).and.
224
     >    (mode.ne.'grid'  ))  
225
     >then
226
         print*,' ERROR: Invalid mode ',trim(mode)
227
         stop
228
      endif
229
 
230
c     Allocate memory for old and new (reparameterised) trajectory
231
      allocate(olon(ntime),stat=stat)
232
      if (stat.ne.0) print*,'*** error allocating array olon ***'
233
      allocate(olat(ntime),stat=stat)
234
      if (stat.ne.0) print*,'*** error allocating array olat ***'
235
      allocate(otim(ntime),stat=stat)
236
      if (stat.ne.0) print*,'*** error allocating array otim ***'
237
      allocate(nlon(1000*ntime),stat=stat)
238
      if (stat.ne.0) print*,'*** error allocating array nlon ***'
239
      allocate(nlat(1000*ntime),stat=stat)
240
      if (stat.ne.0) print*,'*** error allocating array nlat ***'
241
      allocate(ntim(1000*ntime),stat=stat)
242
      if (stat.ne.0) print*,'*** error allocating array ntim ***'
243
      allocate(odist(ntime),stat=stat)
244
      if (stat.ne.0) print*,'*** error allocating array odist ***'
245
      allocate(ndist(1000*ntime),stat=stat)
246
      if (stat.ne.0) print*,'*** error allocating array ndist ***'
247
      allocate(ofld(ntime),stat=stat)
248
      if (stat.ne.0) print*,'*** error allocating array ofld ***'
249
      allocate(nfld(1000*ntime),stat=stat)
250
      if (stat.ne.0) print*,'*** error allocating array nfld ***'
251
 
252
c     Allocate memory for complete trajectory set
253
      allocate(trainp(ntra,ntime,nfield),stat=stat)
254
      if (stat.ne.0) print*,'*** error allocating array trainp ***'
255
      allocate(sel_flag(ntra),stat=stat)
256
      if (stat.ne.0) print*,'*** error allocating array sel_flag ***'
257
 
258
c     Allocate memory for auxiliary fields
259
      allocate(traj(nfield),stat=stat)
260
      if (stat.ne.0) print*,'*** error allocating array traj ***'
261
 
262
c     Set the format of the input file
263
      call mode_tra(inpmode,inpfile)
264
      if (inpmode.eq.-1) inpmode=1
265
 
266
c     Read the input trajectory file
267
      call ropen_tra(fid,inpfile,ntra,ntime,nfield,
268
     >                   reftime,varsinp,inpmode)
269
      call read_tra (fid,trainp,ntra,ntime,nfield,inpmode)
270
      call close_tra(fid,inpmode)
271
 
272
c     Check that first four columns correspond to time,lon,lat,p
273
      if ( (varsinp(1).ne.'time' ).or.
274
     >     (varsinp(2).ne.'xpos' ).and.(varsinp(2).ne.'lon' ).or.
275
     >     (varsinp(3).ne.'ypos' ).and.(varsinp(3).ne.'lat' ).or.
276
     >     (varsinp(4).ne.'ppos' ).and.(varsinp(4).ne.'p'   ) )
277
     >then
278
         print*,' ERROR: problem with input trajectories ...'
279
         stop
280
      endif
281
      varsinp(1) = 'TIME'
282
      varsinp(2) = 'lon'
283
      varsinp(3) = 'lat'
284
      varsinp(4) = 'p'
285
 
286
c     Get the index of the field (if needed)
287
      if ( field.ne.'nil' ) then
288
         ifield = 0
289
         do i=1,nfield
290
            if ( varsinp(i).eq.field ) ifield = i
291
         enddo
292
         if ( ifield.eq.0 ) then
293
            print*,' ERROR: field ',trim(field),' not found... Stop'
294
            stop
295
         endif
296
      endif
297
 
298
c     Write some status information of the input trajectories
299
      print*
300
      print*,'---- INPUT TRAJECTORIES ---------------------------------'
301
      print*
302
      print*,' Reference time (year)  : ',reftime(1)
303
      print*,'                (month) : ',reftime(2)
304
      print*,'                (day)   : ',reftime(3)
305
      print*,'                (hour)  : ',reftime(4)
306
      print*,'                (min)   : ',reftime(5)
307
      print*,' Time range (min)       : ',reftime(6)
308
      do i=1,nfield
309
         if ( i.ne.ifield ) then
310
            print*,' Var                    :',i,trim(varsinp(i))
311
         else
312
            print*,' Var                    :',i,trim(varsinp(i)),
313
     >                                        '       [ gridding ]'
314
         endif
315
      enddo
316
      print*,' List of selected times'
317
      do i=1,ntime
318
         if ( (step.eq.0).or.(step.eq.i) ) then
319
            print*,'     ',i,'  -> ',trainp(1,i,1)
320
         endif
321
      enddo
322
      print*
323
 
324
c     Select flag: all trajectories are selected
325
      if ( sel_file.eq.'nil' ) then
326
 
327
         do i=1,ntra
328
            sel_flag(i) = 1
329
         enddo
330
 
331
c     Select flag: index file
332
      elseif ( sel_format.eq.'index' ) then
333
 
334
         do i=1,ntra
335
            sel_flag(i) = 0
336
         enddo
337
 
338
         open(10,file=sel_file)
339
 142      read(10,*,end=141) ind
340
          sel_flag(ind) = 1
341
          goto 142 
342
 141     continue
343
         close(10)
344
 
345
c     Select flag: boolean file
346
      elseif ( sel_format.eq.'boolean' ) then
347
 
348
         open(10,file=sel_file)
349
          do i=1,ntra
350
            read(10,*) ind
351
            if ( ind.eq.1 ) sel_flag(i) = ind
352
          enddo
353
         close(10)
354
 
355
      endif
356
 
357
c     Write status information
358
      if ( sel_file.eq.'nil' ) then
359
          print*,' Selected trajectories  : all ',ntra          
360
       else
361
          count = 0
362
          do i=1,ntra
363
             if ( sel_flag(i).eq.1 ) count = count + 1
364
          enddo
365
          print*,' #selected trajectories : ',count,
366
     >            ' [ ',real(count)/real(ntra) * 100.,' % ] '
367
       endif
368
       print*
369
 
370
c     ---------------------------------------------------------------------
371
c     Coordinate transformations and grid adjustment
372
c     ---------------------------------------------------------------------
373
 
374
c     Transform from lat/lon to rotated lat/lon, if requested
375
      if ( gridtype.eq.'rotated') then
376
 
377
         crot = 0.
378
 
379
         pollon=clon-180.
380
         if (pollon.lt.-180.) pollon=pollon+360.
381
         pollat=90.-clat
382
         do i=1,ntra
383
            do j=1,ntime
384
 
385
               if ( sel_flag(i).eq.1 ) then
386
 
387
c                Get lat/lon coordinates for trajectory point
388
                 lon = trainp(i,j,2)
389
                 lat = trainp(i,j,3)
390
 
391
c                First Rotation
392
                 pollon=clon-180.
393
                 if (pollon.lt.-180.) pollon=pollon+360.
394
                 pollat=90.-clat
395
                 rlon0=lmtolms(lat,lon,pollat,pollon)
396
                 rlat0=phtophs(lat,lon,pollat,pollon)            
397
 
398
c                Second rotation
399
                 pollon=-180.
400
                 pollat=90.+crot
401
                 rlon=90.+lmtolms(rlat0,rlon0-90.,pollat,pollon)
402
                 rlat=phtophs(rlat0,rlon0-90.,pollat,pollon)   
403
 
404
c                Get rotated latitude and longitude
405
 100             if (rlon.lt.xmin) then
406
                  rlon=rlon+360.
407
                  goto 100
408
                 endif
409
 102             if (rlon.gt.(xmin+real(nx-1)*dx)) then
410
                  rlon=rlon-360.
411
                  goto 102
412
                 endif
413
 
414
c                Set the new trajectory coordinates
415
                 trainp(i,j,2) = rlon
416
                 trainp(i,j,3) = rlat
417
 
418
              endif
419
 
420
            enddo
421
         enddo
422
 
423
      endif
424
 
425
c     Dynamic grid adjustment
426
      if ( dynamic_grid.eq.1 ) then
427
 
428
c        Get the grid parameters
429
         xmin =  180.
430
         ymin =   90.
431
         xmax = -180.
432
         ymax =  -90.
433
 
434
         do i=1,ntra
435
 
436
            if ( sel_flag(i).eq.1 ) then
437
 
438
              if ( step.eq.0 ) then
439
               do j=1,ntime
440
                  if ( trainp(i,j,2).lt.xmin) xmin =  trainp(i,j,2)
441
                  if ( trainp(i,j,2).gt.xmax) xmax =  trainp(i,j,2)
442
                  if ( trainp(i,j,3).lt.ymin) ymin =  trainp(i,j,3)
443
                  if ( trainp(i,j,3).gt.ymax) ymax =  trainp(i,j,3)
444
               enddo
445
              else
446
                if ( trainp(i,step,2).lt.xmin) xmin =  trainp(i,step,2)
447
                if ( trainp(i,step,2).gt.xmax) xmax =  trainp(i,step,2)
448
                if ( trainp(i,step,3).lt.ymin) ymin =  trainp(i,step,3)
449
                if ( trainp(i,step,3).gt.ymax) ymax =  trainp(i,step,3)
450
              endif
451
 
452
            endif
453
 
454
         enddo
455
 
456
c        Get first guess for "optimal" grid
457
         nx = 400
458
         ny = 400
459
         dx = (xmax - xmin)/real(nx-1)
460
         dy = (ymax - ymin)/real(ny-1)
461
 
462
c        Make the grid spacing equal in zonal and meridional direction
463
         if ( dx.gt.dy ) then
464
 
465
            dy = dx
466
            ny = (ymax - ymin)/dy + 1
467
            if (ny.lt.nx/2)              ny = nx / 2
468
            if ( real(ny)*dy .ge. 180. ) ny = 180./dy + 1
469
            ycen = 0.5* (ymin+ymax)
470
            ymin = ycen - 0.5 * real(ny/2) * dy
471
            if (ymin.le.-90.) ymin = -90.
472
 
473
         else
474
 
475
            dx = dy
476
            nx = (xmax - xmin)/dx + 1
477
            if (nx.lt.ny/2)              nx = ny / 2
478
            if ( real(nx)*dx .ge. 360. ) nx = 360./dx + 1
479
            xcen = 0.5* (xmin+xmax)
480
            xmin = xcen - 0.5 * real(nx/2) * dx
481
            if (xmin.le.-180.) xmin = -180.
482
 
483
         endif
484
 
485
c        Write information
486
         print*
487
         print*,'---- DYNAMIC GRID ADJUSTMENT',
488
     >          ' ----------------------------'  
489
         print*
490
         print*,'Grid   nlon,nlat     : ',nx,ny
491
         print*,'       lonmin,latmin : ',xmin,ymin
492
         print*,'       dlon,dlat     : ',dx,dy
493
         print*
494
 
495
c     Write grid information for rotated grid (if not already done
496
      elseif ( gridtype.eq.'rotated') then
497
 
498
         print*
499
         print*,'---- GRID PARAMETERS -------',
500
     >          ' ----------------------------'  
501
         print*
502
         print*,'Grid   nlon,nlat     : ',nx,ny
503
         print*,'       lonmin,latmin : ',xmin,ymin
504
         print*,'       dlon,dlat     : ',dx,dy
505
         print*
506
 
507
 
508
      endif
509
 
510
c     Set the grid boundaries
511
      xmax=xmin+real(nx-1)*dx
512
      ymax=ymin+real(ny-1)*dy
513
 
514
c     Allocate memory for output array and auxiliary gridding array 
515
      allocate(cnt(nx,ny),stat=stat)
516
      if (stat.ne.0) print*,'*** error allocating array cnt  ***'
517
      allocate(res(nx,ny),stat=stat)
518
      if (stat.ne.0) print*,'*** error allocating array res  ***'
519
      allocate(fld(nx,ny),stat=stat)
520
      if (stat.ne.0) print*,'*** error allocating array fld  ***'
521
      allocate(area(nx,ny),stat=stat)
522
      if (stat.ne.0) print*,'*** error allocating array area ***'
523
 
524
      allocate(connect0(nx,ny),stat=stat)
525
      if (stat.ne.0) print*,'*** error allocating array connect0 ***'
526
      allocate(connect1(nx,ny),stat=stat)
527
      if (stat.ne.0) print*,'*** error allocating array connect1 ***'
528
      allocate(connect2(nx,ny),stat=stat)
529
      if (stat.ne.0) print*,'*** error allocating array connect2 ***'
530
 
531
 
532
c     Init the output array
533
      do i=1,nx
534
         do j=1,ny
535
            connect0(i,j) = 0
536
            connect1(i,j) = 0
537
            connect2(i,j) = 0
538
            cnt(i,j)      = 0.
539
            res(i,j)      = 0.
540
            fld(i,j)      = 0.
541
         enddo
542
      enddo  
543
 
544
c     ---------------------------------------------------------------------
545
c     Gridding
546
c     ---------------------------------------------------------------------
547
 
548
c     Write some status information 
549
      print*,'---- GRIDDING -------------------------------------------'
550
      print*
551
 
552
c     Loop over all entries of sampling table
553
      connectval0 = 0
554
      connectval1 = 0
555
      connectval2 = 0
556
      count       = 0
557
 
558
      do i=1,ntra
559
 
560
         if (mod(i,100).eq.0) print*,i,' of ',ntra
561
 
562
c        Skip all trajectories which are not selected
563
         if ( sel_flag(i).eq.0 ) goto 300
564
 
565
c        ------- Read a complete trajectory ---------------------------
566
         do j=1,ntime
567
            otim(j) = trainp(i,j,1)
568
            olon(j) = trainp(i,j,2)
569
            olat(j) = trainp(i,j,3)
570
            if ( field.ne.'nil' ) then
571
               ofld(j) =trainp(i,j,ifield)
572
            endif
573
         enddo
574
 
575
c        -------- Convert hh.m time into fractional time --------------
576
         do j=1,ntime
577
            hhmm    = otim(j)
578
            call hhmm2frac (hhmm,frac)
579
            otim(j) = frac
580
         enddo
581
 
582
c        -------- Interpolation ---------------------------------------
583
 
584
c        Keep the trajectory points as they are
585
         if ( ( mode.eq.'keep').and.(step.eq.0) ) then
586
            npts=opts
587
            do j=1,opts
588
               ntim(j)=otim(j)
589
               nlon(j)=olon(j)
590
               nlat(j)=olat(j)
591
               if ( field.ne.'nil' ) then
592
                  nfld(j)=ofld(j)
593
               endif
594
            enddo
595
 
596
c        Select a single time step
597
         elseif ( ( mode.eq.'keep').and.(step.gt.0) ) then
598
            npts    = 1
599
            ntim(1) = otim(step)
600
            nlon(1) = olon(step)
601
            nlat(1) = olat(step)
602
            if ( field.ne.'nil' ) then
603
               nfld(1) = ofld(step)
604
            endif
605
 
606
c        Perform a reparameterisation in time
607
         else if ( (mode.eq.'time').and.(step.eq.0) ) then
608
 
609
c           Get the new number of trajectory points
610
            npts=nint(abs(otim(opts)-otim(1))/param)+1
611
 
612
c           Handle date line problem
613
            do j=2,opts
614
               if ( (olon(j-1)-olon(j)).gt.180. ) then
615
                  olon(j) = olon(j) + 360.
616
               else if ( (olon(j-1)-olon(j)).lt.-180. ) then
617
                  olon(j) = olon(j) - 360.
618
               endif
619
            enddo
620
 
621
c           Cubic spline fitting
622
            call curvefit(otim,olon,opts,ntim,nlon,npts)
623
            call curvefit(otim,olat,opts,ntim,nlat,npts)
624
            if ( field.ne.'nil' ) then
625
               call curvefit(otim,ofld,opts,ntim,nfld,npts)
626
            endif
627
 
628
c           Reverse date line handling
629
            do j=1,npts
630
               if ( nlon(j).gt.xmax ) then
631
                  nlon(j) = nlon(j) -360.
632
               else if ( nlon(j).lt.xmin ) then
633
                  nlon(j) = nlon(j) +360.
634
               endif
635
            enddo
636
 
637
c        Perform a reparameterisation with equally spaced gridpoint
638
         elseif ( (mode.eq.'space').and.(step.eq.0) ) then
639
 
640
c           Calculate the distance and spacing
641
            odist(1) = 0.
642
            unit     = 'km'
643
            do j=2,ntime
644
               odist(j)=odist(j-1) + 
645
     >                  sdis(olon(j-1),olat(j-1),olon(j),olat(j),unit)
646
            enddo
647
 
648
c           Determine the new number of trajectory points
649
            npts=nint(odist(ntime)/param)+1
650
            if (npts.eq.0) then
651
               npts=1.
652
            endif
653
 
654
c           Handle date line problem
655
            do j=2,opts
656
               if ( (olon(j-1)-olon(j)).gt.180. ) then
657
                  olon(j) = olon(j) + 360.
658
               else if ( (olon(j-1)-olon(j)).lt.-180. ) then
659
                  olon(j) = olon(j) - 360.
660
               endif
661
            enddo
662
 
663
c           Cubic spline fitting
664
            call curvefit(odist,olon,opts,ndist,nlon,npts)
665
            call curvefit(odist,olat,opts,ndist,nlat,npts)
666
            call curvefit(odist,otim,opts,ndist,ntim,npts)
667
            if ( field.ne.'nil' ) then
668
               call curvefit(odist,ofld,opts,ndist,nfld,npts)
669
            endif
670
 
671
c           Reverse date line handling
672
            do j=1,npts
673
               if ( nlon(j).gt.xmax ) then
674
                  nlon(j) = nlon(j) -360.
675
               else if ( nlon(j).lt.xmin ) then
676
                  nlon(j) = nlon(j) +360.
677
               endif
678
            enddo
679
 
680
c        Perform a reparameterisation with equally spaced gridpoint
681
         elseif ( (mode.eq.'grid').and.(step.eq.0) ) then
682
 
683
c           Calculate the distance and spacing
684
            odist(1) = 0.
685
            unit     = 'deg'
686
            do j=2,ntime
687
               odist(j)=odist(j-1) + 
688
     >                  sdis(olon(j-1),olat(j-1),olon(j),olat(j),unit)
689
            enddo
690
 
691
c           Determine the new number of trajectory points
692
            npts=nint(odist(ntime)/param)+1
693
            if (npts.eq.0) then
694
               npts=1.
695
            endif
696
 
697
c           Handle date line problem
698
            do j=2,opts
699
               if ( (olon(j-1)-olon(j)).gt.180. ) then
700
                  olon(j) = olon(j) + 360.
701
               else if ( (olon(j-1)-olon(j)).lt.-180. ) then
702
                  olon(j) = olon(j) - 360.
703
               endif
704
            enddo
705
 
706
c           Cubic spline fitting
707
            call curvefit(odist,olon,opts,ndist,nlon,npts)
708
            call curvefit(odist,olat,opts,ndist,nlat,npts)
709
            call curvefit(odist,otim,opts,ndist,ntim,npts)
710
            if ( field.ne.'nil' ) then
711
               call curvefit(odist,ofld,opts,ndist,nfld,npts)
712
            endif
713
 
714
c           Reverse date line handling
715
            do j=1,npts
716
               if ( nlon(j).gt.xmax ) then
717
                  nlon(j) = nlon(j) -360.
718
               else if ( nlon(j).lt.xmin ) then
719
                  nlon(j) = nlon(j) +360.
720
               endif
721
            enddo
722
 
723
         endif
724
 
725
c        -------- Do the gridding -------------------------------------
726
 
727
c        Gridding of trajectory
728
         do j=1,npts
729
 
730
c           Check whether point is in data domain
731
	    if ( (nlon(j).gt.xmin).and.(nlon(j).lt.xmax).and.
732
     >           (nlat(j).gt.ymin).and.(nlat(j).lt.ymax))
733
     >      then
734
 
735
c              Increase counter for gridded points
736
               count = count + 1
737
 
738
c              ----------------- Gridding: simple count -----------------
739
               connectval0 = connectval0+1
740
 
741
               addvalue    = 1.
742
 
743
               call  gridding1
744
     >              (nlat(j),nlon(j),addvalue,
745
     >               radius,runit,connect0,connectval0,
746
     >               cnt,nx,ny,xmin,ymin,dx,dy)
747
 
748
c              ----------------- Gridding: residence time ---------------
749
               connectval1 = connectval1+1
750
 
751
               if ( ntime.eq.1 ) then
752
                  addvalue = 0.
753
               elseif ( j.eq.1 )  then
754
                  addvalue=abs(ntim(2)-ntim(1))
755
               else
756
                  addvalue=abs(ntim(j)-ntim(j-1))
757
               endif
758
 
759
               call  gridding1
760
     >              (nlat(j),nlon(j),addvalue,
761
     >               radius,runit,connect1,connectval1,
762
     >               res,nx,ny,xmin,ymin,dx,dy)
763
 
764
 
765
c              --------------- Gridding: field -------------------------
766
               if ( field.ne.'nil' ) then
767
 
768
                   connectval2 = connectval2+1
769
 
770
                   addvalue    = nfld(j)
771
 
772
                   call  gridding1
773
     >                  (nlat(j),nlon(j),addvalue,
774
     >                  radius,runit,connect2,connectval2,
775
     >                  fld,nx,ny,xmin,ymin,dx,dy)
776
 
777
               endif
778
 
779
	    endif
780
 
781
         enddo
782
 
783
c        Exit point for loop over all trajectories
784
 300     continue
785
 
786
      enddo
787
 
788
c     Write status information
789
      print*
790
      print*,' # gridded points       : ',count
791
 
792
c     ---------------------------------------------------------------------
793
c     Unit conversions and output to netCDF file
794
c     ---------------------------------------------------------------------
795
 
796
c     Write some status information 
797
      print*
798
      print*,'---- WRITE OUTPUT ---------------------------------------'
799
      print*
800
 
801
c     Area (in km^2)
802
      do i=1,nx	         
803
         do j=1,ny	
804
            slat=ymin+real(j-1)*dy
805
            if (abs(abs(slat)-90.).gt.eps) then
806
               area(i,j) = dy*dx*cos(pi180*slat)*deltay**2
807
            else
808
               area(i,j) = 0.
809
            endif
810
         enddo
811
      enddo
812
 
813
c     Normalise gridded field
814
      if ( field.ne.'nil' ) then
815
         do i=1,nx
816
            do j=1,ny
817
               if ( cnt(i,j).gt.0. ) then
818
                  fld(i,j) = fld(i,j) / cnt(i,j)
819
               endif
820
            enddo
821
         enddo
822
      endif
823
 
824
c     Set the time for the output netCDF files - if a composite is
825
c     calculatd, then the time is set to 
826
      if ( step.eq.0 ) then
827
         time = -999.
828
         print*,'   ... COMPOSITE OVER ALL TRAJECTORY TIMES (-999)'
829
         print*
830
      else
831
         time = trainp(1,step,1)
832
      endif
833
 
834
c     Write output to CF netCDF
835
      cdfname  = outfile
836
 
837
      varname  = 'COUNT'
838
      longname = 'trajectory counts'
839
      varunit  = 'counts per grid point'
840
      call  writecdf2D_cf (cdfname,varname,longname,varunit,gridtype,
841
     >       clon,clat,nlonlat,dlonlat,cnt,time,dx,dy,xmin,ymin,nx,
842
     >       ny,crefile,crefile,1)
843
      write(*,'(a8,a10,a5,a10,a10,f7.2,a2)') 
844
     >     '    ... ',trim(varname),' -> ',trim(cdfname),
845
     >     ' [ time = ',time,' ]'    
846
 
847
      varname  = 'RESIDENCE'
848
      longname = 'residence time'
849
      varunit  = 'hours per grid point'
850
      call  writecdf2D_cf (cdfname,varname,longname,varunit,gridtype,
851
     >       clon,clat,nlonlat,dlonlat,res,time,dx,dy,xmin,ymin,nx,
852
     >       ny,0,crefile,1)
853
      write(*,'(a8,a10,a5,a10,a10,f7.2,a2)') 
854
     >     '    ... ',trim(varname),' -> ',trim(cdfname),
855
     >     ' [ time = ',time,' ]'    
856
 
857
      varname  = 'AREA'
858
      longname = 'area corresponding to grid points'
859
      varunit  = 'square kilometers'
860
      call  writecdf2D_cf (cdfname,varname,longname,varunit,gridtype,
861
     >       clon,clat,nlonlat,dlonlat,area,time,dx,dy,xmin,ymin,nx,
862
     >       ny,0,crefile,1)
863
      write(*,'(a8,a10,a5,a10,a10,f7.2,a2)') 
864
     >     '    ... ',trim(varname),' -> ',trim(cdfname),
865
     >     ' [ time = ',time,' ]'    
866
 
867
      if ( field.ne.'nil' ) then
868
         varname  = field
869
         longname = field
870
         varunit  = 'as on trajectory file'
871
         call  writecdf2D_cf (cdfname,varname,longname,varunit,gridtype,
872
     >       clon,clat,nlonlat,dlonlat,fld,time,dx,dy,xmin,ymin,nx,
873
     >       ny,0,crevar,1)
874
 
875
         write(*,'(a8,a10,a5,a10,a10,f7.2,a2)') 
876
     >        '    ... ',trim(varname),' -> ',trim(cdfname),
877
     >        ' [ time = ',time,' ]'    
878
      endif
879
 
880
c     Write status information
881
      print*
882
      print*,'              *** END OF PROGRAM DENSITY **'
883
      print*,'========================================================='
884
 
885
      end
886
 
887
c     ********************************************************************
888
c     * GRIDDING SUBROUTINES                                             *
889
c     ********************************************************************
890
 
891
c     ---------------------------------------------------------------------
892
c     Gridding of one single data point (smoothing in km, deg, gridp)
893
c     ---------------------------------------------------------------------
894
 
895
      subroutine gridding1 (lat,lon,addval,radius,unit,
896
     >                      connect,connectval,
897
     >                      out,nx,ny,xmin,ymin,dx,dy)
898
 
899
      implicit none
900
 
901
c     Declaration of subroutine parameters
902
      real         lat,lon
903
      integer      nx,ny
904
      real         xmin,ymin,dx,dy
905
      real         out(nx,ny)
906
      real         radius
907
      character*80 unit
908
      integer      connectval
909
      integer      connect(nx,ny)
910
      real         addval
911
 
912
c     Auxiliary variables
913
      integer   i,j,k
914
      integer   mu,md,nr,nl,n,m
915
      integer   stackx(nx*ny),stacky(nx*ny)
916
      integer   tab_x(nx*ny),tab_y(nx*ny)
917
      real      tab_r(nx*ny)
918
      integer   sp
919
      real      lat2,lon2
920
      real      dist,sum
921
      real      xmax
922
      integer   periodic
923
      integer   test
924
 
925
c     Numerical epsilon
926
      real      eps
927
      parameter (eps=0.01)
928
 
929
c     Externals
930
      real      sdis,weight
931
      external  sdis,weight
932
 
933
c     Check whether lat/lon point is valid
934
      xmax=xmin+real(nx-1)*dx
935
      if (lon.lt.xmin-eps) lon=lon+360.
936
      if (lon.gt.xmax+eps) lon=lon-360.
937
      if (abs(lat-90).lt.eps) lat=90.
938
      if (abs(lat+90).lt.eps) lat=-90.
939
      if ((abs(lat).gt.(90.+eps)).or.
940
     >    (lon.lt.xmin-eps).or.(lon.gt.xmax+eps)) then
941
         print*,'Invalid lat/lon point ',lat,lon
942
         return
943
      endif
944
 
945
c     Set flag for periodic domain
946
      if (abs(xmax-xmin-360.).lt.eps) then
947
         periodic=1
948
      else if (abs(xmax-xmin-360+dx).lt.eps) then
949
         periodic=2
950
      else
951
         periodic=0
952
      endif
953
 
954
c     Get indices of one coarse grid point within search radius
955
      i=nint((lon-xmin)/dx)+1
956
      if ((i.eq.nx).and.(periodic.eq.1)) i=1
957
      j=nint((lat-ymin)/dy)+1
958
      lat2=ymin+real(j-1)*dy
959
      lon2=xmin+real(i-1)*dx
960
      dist=sdis(lon,lat,lon2,lat2,unit)
961
      if (dist.gt.radius) then
962
         print*,'1: Search radius is too small...'
963
         stop
964
      endif
965
 
966
c     Get connected points
967
      k=0
968
      stackx(1)=i
969
      stacky(1)=j
970
      sp=1
971
      do while (sp.ne.0) 
972
 
973
c        Get an element from stack
974
         n=stackx(sp)
975
         m=stacky(sp)
976
         sp=sp-1
977
 
978
c        Get distance from reference point
979
         lat2=ymin+real(m-1)*dy
980
         lon2=xmin+real(n-1)*dx
981
         dist=sdis(lon,lat,lon2,lat2,unit)
982
 
983
c        Check whether distance is smaller than search radius: connected
984
         if (dist.lt.radius) then
985
 
986
c           Make entry in filter mask
987
            k=k+1
988
            tab_x(k)=n
989
            tab_y(k)=m
990
            tab_r(k)=weight(dist,radius)
991
 
992
c           Mark this point as visited
993
            connect(n,m)=connectval
994
 
995
c           Get coordinates of neighbouring points
996
            nr=n+1
997
            if ((nr.gt.nx)  .and.(periodic.eq.0)) nr=nx
998
            if ((nr.gt.nx-1).and.(periodic.eq.1)) nr=1
999
            if ((nr.gt.nx)  .and.(periodic.eq.2)) nr=1
1000
            nl=n-1
1001
            if ((nl.lt.1).and.(periodic.eq.0)) nl=1
1002
            if ((nl.lt.1).and.(periodic.eq.1)) nl=nx-1
1003
            if ((nl.lt.1).and.(periodic.eq.2)) nl=nx
1004
            mu=m+1
1005
            if (mu.gt.ny) mu=ny
1006
            md=m-1
1007
            if (md.lt.1) md=1
1008
 
1009
c           Update stack
1010
            if (connect(nr,m).ne.connectval) then
1011
               connect(nr,m)=connectval
1012
               sp=sp+1
1013
               stackx(sp)=nr
1014
               stacky(sp)=m
1015
            endif
1016
            if (connect(nl,m).ne.connectval) then
1017
               connect(nl,m)=connectval
1018
               sp=sp+1
1019
               stackx(sp)=nl
1020
               stacky(sp)=m
1021
            endif
1022
            if (connect(n,mu).ne.connectval) then
1023
               connect(n,mu)=connectval
1024
               sp=sp+1
1025
               stackx(sp)=n
1026
               stacky(sp)=mu
1027
            endif
1028
            if (connect(n,md).ne.connectval) then
1029
               connect(n,md)=connectval
1030
               sp=sp+1
1031
               stackx(sp)=n
1032
               stacky(sp)=md
1033
            endif
1034
         endif
1035
 
1036
      end do
1037
 
1038
      if (k.ge.1) then
1039
         sum=0.
1040
         do i=1,k
1041
            sum=sum+tab_r(i)
1042
         enddo
1043
         do i=1,k
1044
            out(tab_x(i),tab_y(i))=out(tab_x(i),tab_y(i))+
1045
     >                             addval*tab_r(i)/sum
1046
 
1047
            if ((tab_x(i).eq.1).and.(periodic.eq.1)) then
1048
               out(nx,tab_y(i))=out(nx,tab_y(i))+
1049
     >                             addval*tab_r(i)/sum
1050
            endif
1051
         enddo
1052
      else
1053
         print*,'2: Search radius is too small...'
1054
         stop
1055
      endif
1056
 
1057
      end
1058
 
1059
 
1060
c     ----------------------------------------------------------------------
1061
c     Get spherical distance between lat/lon points
1062
c     ----------------------------------------------------------------------
1063
 
1064
      real function sdis(xp,yp,xq,yq,unit)
1065
 
1066
c     Calculates spherical distance (in km) between two points given
1067
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
1068
 
1069
      real         re
1070
      parameter    (re=6370.)
1071
      real         xp,yp,xq,yq,arg
1072
      character*80 unit
1073
      real         dlon
1074
 
1075
      if ( unit.eq.'km' ) then
1076
 
1077
         arg=sind(yp)*sind(yq)+cosd(yp)*cosd(yq)*cosd(xp-xq)
1078
         if (arg.lt.-1.) arg=-1.
1079
         if (arg.gt.1.) arg=1.
1080
         sdis=re*acos(arg)
1081
 
1082
      elseif ( unit.eq.'deg' ) then
1083
 
1084
         dlon = xp-xq
1085
         if ( dlon.gt. 180. ) dlon = dlon - 360.
1086
         if ( dlon.lt.-180. ) dlon = dlon + 360.
1087
         sdis = sqrt( dlon**2 + (yp-yq)**2 )
1088
 
1089
      endif
1090
 
1091
 
1092
c     Quick and dirty trick to avoid zero distances
1093
      if (sdis.eq.0.) sdis=0.1
1094
 
1095
      end
1096
 
1097
c     ----------------------------------------------------------------------
1098
c     Weight function for the filter mask
1099
c     ----------------------------------------------------------------------
1100
 
1101
      real function weight (r,radius)
1102
 
1103
c     Attribute to each distanc r its corresponding weight in the filter mask
1104
 
1105
      implicit none
1106
 
1107
c     Declaration of subroutine parameters
1108
      real r
1109
      real radius
1110
 
1111
c     Simple 0/1 mask
1112
      if (r.lt.radius) then
1113
         weight=exp(-r/radius)
1114
      else
1115
         weight=0.
1116
      endif
1117
 
1118
      end
1119
 
1120
 
1121
c     ********************************************************************
1122
c     * REPARAMETERIZATION SUBROUTINES                                   *
1123
c     ********************************************************************
1124
 
1125
c     -------------------------------------------------------------
1126
c     Interpolation of the trajectory with a natural cubic spline
1127
c     -------------------------------------------------------------
1128
 
1129
      SUBROUTINE curvefit (time,lon,n,
1130
     >                     sptime,splon,spn)
1131
 
1132
c     Given the curve <time,lon> with <n> data points, fit a
1133
c     cubic spline to this curve. The new curve is returned in 
1134
c     <sptime,splon,spn> with <spn> data points. The parameter
1135
c     <spn> specifies on entry the number of spline interpolated points
1136
c     along the curve.
1137
 
1138
      implicit none
1139
 
1140
c     Declaration of subroutine parameters
1141
      integer n
1142
      real time(n),lon(n)
1143
      integer spn
1144
      real sptime(spn),splon(spn)
1145
 
1146
c     Auxiliary variables
1147
      real y2ax(n)
1148
      real dt
1149
      real s
1150
      integer i
1151
      real order
1152
 
1153
c     Determine whether the input array is ascending or descending
1154
      if (time(1).gt.time(n)) then
1155
         order=-1.
1156
      else
1157
         order= 1.
1158
      endif
1159
 
1160
c     Bring the time array into ascending order
1161
      do i=1,n
1162
         time(i)=order*time(i)
1163
      enddo
1164
 
1165
c     Prepare the (natural) cubic spline interpolation
1166
      call spline (time,lon,n,1.e30,1.e30,y2ax)
1167
      dt=(time(n)-time(1))/real(spn-1)
1168
      do i=1,spn
1169
         sptime(i)=time(1)+real(i-1)*dt
1170
      enddo
1171
 
1172
c     Do the spline interpolation
1173
      do i=1,spn
1174
         call splint(time,lon,y2ax,n,sptime(i),s)
1175
         splon(i)=s
1176
      enddo
1177
 
1178
c     Change the time arrays back
1179
      do i=1,spn
1180
         sptime(i)=order*sptime(i)
1181
      enddo
1182
      do i=1,n
1183
         time(i)=order*time(i)
1184
      enddo
1185
 
1186
      return
1187
      end
1188
 
1189
c     -------------------------------------------------------------
1190
c     Basic routines for spline interpolation (Numerical Recipes)
1191
c     -------------------------------------------------------------
1192
 
1193
      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
1194
      INTEGER n,NMAX
1195
      REAL yp1,ypn,x(n),y(n),y2(n)
1196
      PARAMETER (NMAX=500)
1197
      INTEGER i,k
1198
      REAL p,qn,sig,un,u(NMAX)
1199
      if (yp1.gt..99e30) then
1200
        y2(1)=0.
1201
        u(1)=0.
1202
      else
1203
        y2(1)=-0.5
1204
        u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
1205
      endif
1206
      do 11 i=2,n-1
1207
        sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
1208
        p=sig*y2(i-1)+2.
1209
        y2(i)=(sig-1.)/p
1210
        u(i)=(6.*((y(i+1)-y(i))/(x(i+
1211
     *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*
1212
     *u(i-1))/p
1213
11    continue
1214
      if (ypn.gt..99e30) then
1215
        qn=0.
1216
        un=0.
1217
      else
1218
        qn=0.5
1219
        un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
1220
      endif
1221
      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.)
1222
      do 12 k=n-1,1,-1
1223
        y2(k)=y2(k)*y2(k+1)+u(k)
1224
12    continue
1225
      return
1226
      END
1227
 
1228
      SUBROUTINE splint(xa,ya,y2a,n,x,y)
1229
      INTEGER n
1230
      REAL x,y,xa(n),y2a(n),ya(n)
1231
      INTEGER k,khi,klo
1232
      REAL a,b,h
1233
      klo=1
1234
      khi=n
1235
1     if (khi-klo.gt.1) then
1236
        k=(khi+klo)/2
1237
        if(xa(k).gt.x)then
1238
          khi=k
1239
        else
1240
          klo=k
1241
        endif
1242
      goto 1
1243
      endif
1244
      h=xa(khi)-xa(klo)
1245
      if (h.eq.0.) pause 'bad xa input in splint'
1246
      a=(xa(khi)-x)/h
1247
      b=(x-xa(klo))/h
1248
      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h**
1249
     *2)/6.
1250
      return
1251
      END
1252
 
1253
c     ********************************************************************
1254
c     * INPUT / OUTPUT SUBROUTINES                                       *
1255
c     ********************************************************************
1256
 
1257
 
1258
c     --------------------------------------------------------------------
1259
c     Subroutines to write the CF netcdf output file
1260
c     --------------------------------------------------------------------
1261
 
1262
      subroutine writecdf2D_cf 
1263
     >          (cdfname,varname,longname,unit,gridtype,clon,clat,
1264
     >           nlonlat,dlonlat,arr,time,dx,dy,xmin,ymin,nx,ny,
1265
     >           crefile,crevar,cretime)
1266
 
1267
c     Create and write to the CF netcdf file <cdfname>. The variable
1268
c     with name <varname> and with time <time> is written. The data
1269
c     are in the two-dimensional array <arr>. The list <dx,dy,xmin,
1270
c     ymin,nx,ny> specifies the output grid. The flags <crefile> and
1271
c     <crevar> determine whether the file and/or the variable should
1272
c     be created; correspondingly for the unlimited dimension <time>
1273
c     with the flag <cretime>.
1274
 
1275
      USE netcdf
1276
 
1277
      IMPLICIT NONE
1278
 
1279
c     Declaration of input parameters
1280
      character*80 cdfname
1281
      character*80 varname,longname,unit
1282
      integer      nx,ny
1283
      real         arr(nx,ny)
1284
      real         dx,dy,xmin,ymin
1285
      real         time
1286
      integer      crefile,crevar,cretime
1287
      character*80 gridtype
1288
      real         clon,clat
1289
      integer      nlonlat
1290
      real         dlonlat
1291
 
1292
c     Local variables
1293
      integer      ierr
1294
      integer      ncID
1295
      integer      LonDimId,    varLonID
1296
      integer      LatDimID,    varLatID
1297
      integer      TimeDimID,   varTimeID
1298
      real         longitude(nx)
1299
      real         latitude (ny)
1300
      real         timeindex
1301
      integer      i
1302
      integer      nvars,varids(100)
1303
      integer      ndims,dimids(100)
1304
      real         timelist(1000)
1305
      integer      ntimes
1306
      integer      ind
1307
      integer      varID
1308
 
1309
c     Quick an dirty solution for fieldname conflict
1310
      if ( varname.eq.'time' ) varname = 'TIME'
1311
 
1312
c     Initially set error to indicate no errors.
1313
      ierr = 0
1314
 
1315
c     ---- Create the netCDF - skip if <crefile=0> ----------------------
1316
      if ( crefile.ne.1 ) goto 100
1317
 
1318
c     Create the file 
1319
      ierr = nf90_create(trim(cdfname), NF90_CLOBBER, ncID)
1320
 
1321
c     Define dimensions 
1322
      ierr=nf90_def_dim(ncID,'longitude',nx            , LonDimID )
1323
      ierr=nf90_def_dim(ncID,'latitude' ,ny            , LatDimID )
1324
      ierr=nf90_def_dim(ncID,'time'     ,nf90_unlimited, TimeDimID)
1325
 
1326
c     Define coordinate Variables 
1327
      ierr = nf90_def_var(ncID,'longitude',NF90_FLOAT,
1328
     >     (/ LonDimID /),varLonID)
1329
      ierr = nf90_put_att(ncID, varLonID, "standard_name","longitude")
1330
      ierr = nf90_put_att(ncID, varLonID, "units"      ,"degree_east")
1331
 
1332
      ierr = nf90_def_var(ncID,'latitude',NF90_FLOAT,
1333
     >     (/ LatDimID /),varLatID)
1334
      ierr = nf90_put_att(ncID, varLatID, "standard_name", "latitude")
1335
      ierr = nf90_put_att(ncID, varLatID, "units"    ,"degree_north")
1336
 
1337
      ierr = nf90_def_var(ncID,'time',NF90_FLOAT, 
1338
     >     (/ TimeDimID /), varTimeID)
1339
      ierr = nf90_put_att(ncID, varTimeID, "axis",            "T")
1340
      ierr = nf90_put_att(ncID, varTimeID, "calendar", "standard")
1341
      ierr = nf90_put_att(ncID, varTimeID, "long_name",    "time")
1342
      ierr = nf90_put_att(ncID, varTimeID, "units",       "hours")
1343
 
1344
c     Write global attributes 
1345
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'Conventions', 'CF-1.0')
1346
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'title',  
1347
     >     'Trajectory Densities')
1348
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'source', 
1349
     >     'Lagranto Trajectories')
1350
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'institution', 
1351
     >     'ETH Zurich, IACETH')
1352
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'grid',trim(gridtype) ) 
1353
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'clon',clon )
1354
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'clat',clat )
1355
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'nlonlat',nlonlat )
1356
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'dlonlat',dlonlat )
1357
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'nx',nx )
1358
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ny',ny )
1359
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'dx',dx )
1360
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'dy',dy )
1361
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'xmin',xmin )
1362
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'ymin',ymin )
1363
 
1364
c     Write coordinate data 
1365
      do i = 1,nx+1
1366
         longitude(i) = xmin + real(i-1) * dx
1367
      enddo
1368
      do i = 1,ny+1
1369
         latitude(i)  = ymin + real(i-1) * dy
1370
      enddo
1371
 
1372
c     Check whether the definition was successful
1373
      ierr = nf90_enddef(ncID)
1374
      if (ierr.gt.0) then
1375
         print*, 'An error occurred while attempting to ', 
1376
     >        'finish definition mode.'
1377
         stop
1378
      endif
1379
 
1380
c     Write coordinate data  
1381
      ierr = nf90_put_var(ncID,varLonID ,longitude)
1382
      ierr = nf90_put_var(ncID,varLatID ,latitude )
1383
 
1384
c     Close netCDF file 
1385
      ierr = nf90_close(ncID)
1386
 
1387
 100  continue
1388
 
1389
c     ---- Define a new variable - skip if <crevar=0> -----------------------
1390
 
1391
      if ( crevar.ne.1 ) goto 110
1392
 
1393
c     Open the file for read(write access
1394
      ierr = nf90_open  (trim(cdfname), NF90_WRITE  , ncID)
1395
 
1396
c     Get the IDs for dimensions
1397
      ierr = nf90_inq_dimid(ncID,'longitude', LonDimID )
1398
      ierr = nf90_inq_dimid(ncID,'latitude' , LatDimID )
1399
      ierr = nf90_inq_dimid(ncID,'time'     , TimeDimID)
1400
 
1401
c     Enter define mode
1402
      ierr = nf90_redef(ncID)
1403
 
1404
c     Write definition and add attributes
1405
      ierr = nf90_def_var(ncID,varname,NF90_FLOAT,
1406
     >                   (/ LonDimID, LatDimID, varTimeID /),varID)
1407
      ierr = nf90_put_att(ncID, varID, "long_name" , longname )
1408
      ierr = nf90_put_att(ncID, varID, "units"     , unit     ) 
1409
      ierr = nf90_put_att(ncID, varID, '_FillValue', -999.99  ) 
1410
 
1411
c     Check whether definition was successful
1412
      ierr = nf90_enddef(ncID)
1413
      if (ierr.gt.0) then
1414
         print*, 'An error occurred while attempting to ', 
1415
     >           'finish definition mode.'
1416
         stop
1417
      endif
1418
 
1419
c     Close netCDF file 
1420
      ierr = nf90_close(ncID)
1421
 
1422
 110  continue
1423
 
1424
c     ---- Create a new time (unlimited dimension) - skip if <cretime=0> ------
1425
 
1426
      if ( cretime.ne.1 ) goto 120
1427
 
1428
c     Open the file for read/write access
1429
      ierr = nf90_open  (trim(cdfname), NF90_WRITE, ncID)
1430
 
1431
c     Get the list of times on the netCDF file
1432
      ierr = nf90_inq_dimid(ncID,'time', TimeDimID)
1433
      if ( ierr.ne.0 ) then
1434
         print*,'Time dimension is not defined on ',trim(cdfname),
1435
     >          ' .... Stop'
1436
         stop
1437
      endif
1438
      ierr = nf90_inquire_dimension(ncID, TimeDimID, len = ntimes)
1439
      ierr = nf90_inq_varid(ncID,'time', varTimeID)
1440
      if ( ierr.ne.0 ) then
1441
         print*,'Variable time is not defined on ',trim(cdfname),
1442
     >          ' ... Stop'
1443
         stop
1444
      endif
1445
      ierr = nf90_get_var(ncID,varTimeID,timelist(1:ntimes))
1446
 
1447
c     Decide whether a new time must be written
1448
      ind = 0
1449
      do i=1,ntimes
1450
         if ( time.eq.timelist(i) ) ind = i
1451
      enddo
1452
 
1453
c     Extend the time list if required 
1454
      if ( ind.eq.0 ) then
1455
         ntimes           = ntimes + 1
1456
         timelist(ntimes) = time
1457
         ierr = nf90_put_var(ncID,varTimeID,timelist(1:ntimes))
1458
      endif
1459
 
1460
c     Close netCDF file 
1461
      ierr = nf90_close(ncID)
1462
 
1463
 120  continue
1464
 
1465
c     ---- Write data --------------------------------------------------
1466
 
1467
c     Open the file for read/write access
1468
      ierr = nf90_open  (trim(cdfname), NF90_WRITE , ncID)
1469
 
1470
c     Get the varID
1471
      ierr = nf90_inq_varid(ncID,varname, varID )
1472
      if (ierr.ne.0) then
1473
         print*,'Variable ',trim(varname),' is not defined on ',
1474
     >          trim(cdfname)
1475
         stop
1476
      endif
1477
 
1478
c     Get the time index
1479
      ierr = nf90_inq_dimid(ncID,'time', TimeDimID)
1480
      if ( ierr.ne.0 ) then
1481
         print*,'Time dimension is not defined on ',trim(cdfname),
1482
     >          ' .... Stop'
1483
         stop
1484
      endif
1485
      ierr = nf90_inquire_dimension(ncID, TimeDimID, len = ntimes)
1486
      ierr = nf90_inq_varid(ncID,'time', varTimeID)
1487
      if ( ierr.ne.0 ) then
1488
         print*,'Variable time is not defined on ',trim(cdfname),
1489
     >          ' ... Stop'
1490
         stop
1491
      endif
1492
      ierr = nf90_get_var(ncID,varTimeID,timelist(1:ntimes))
1493
      ind = 0
1494
      do i=1,ntimes
1495
         if ( time.eq.timelist(i) ) ind = i
1496
      enddo
1497
      if (ind.eq.0) then
1498
         print*,'Time',time,' is not defined on the netCDF file',
1499
     >          trim(cdfname),' ... Stop'
1500
         stop
1501
      endif
1502
 
1503
c     Write data block      
1504
      ierr = nf90_put_var(ncID,varID,arr,
1505
     >                    start = (/ 1, 1, ind /), 
1506
     >                    count = (/ nx, ny, 1 /) )
1507
 
1508
c     Check whether writing was successful 
1509
      ierr = nf90_close(ncID)
1510
      if (ierr.ne.0) then
1511
         write(*,*) trim(nf90_strerror(ierr))
1512
         write(*,*) 'An error occurred while attempting to ', 
1513
     >              'close the netcdf file.'
1514
         write(*,*) 'in clscdf_CF'
1515
      endif
1516
 
1517
      end
1518
 
1519
 
1520
c     ********************************************************************************
1521
c     * Transformation routine: LMSTOLM and PHSTOPH from library gm2em               *
1522
c     ********************************************************************************
1523
 
1524
      REAL FUNCTION LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
1525
C
1526
C**** LMSTOLM  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
1527
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
1528
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
1529
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1530
C**   AUFRUF   :   LAM = LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
1531
C**   ENTRIES  :   KEINE
1532
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
1533
C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
1534
C**                IM ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1535
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1536
C**   VERSIONS-
1537
C**   DATUM    :   03.05.90
1538
C**
1539
C**   EXTERNALS:   KEINE
1540
C**   EINGABE-
1541
C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
1542
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
1543
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
1544
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
1545
C**   AUSGABE-
1546
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
1547
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
1548
C**
1549
C**   COMMON-
1550
C**   BLOECKE  :   KEINE
1551
C**
1552
C**   FEHLERBE-
1553
C**   HANDLUNG :   KEINE
1554
C**   VERFASSER:   D.MAJEWSKI
1555
 
1556
      REAL        LAMS,PHIS,POLPHI,POLLAM
1557
 
1558
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
1559
 
1560
      ZSINPOL = SIN(ZPIR18*POLPHI)
1561
      ZCOSPOL = COS(ZPIR18*POLPHI)
1562
      ZLAMPOL = ZPIR18*POLLAM
1563
      ZPHIS   = ZPIR18*PHIS
1564
      ZLAMS   = LAMS
1565
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
1566
      ZLAMS   = ZPIR18*ZLAMS
1567
 
1568
      ZARG1   = SIN(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
1569
     1                          ZCOSPOL*           SIN(ZPHIS)) -
1570
     2          COS(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
1571
      ZARG2   = COS(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
1572
     1                          ZCOSPOL*           SIN(ZPHIS)) +
1573
     2          SIN(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
1574
      IF (ABS(ZARG2).LT.1.E-30) THEN
1575
        IF (ABS(ZARG1).LT.1.E-30) THEN
1576
          LMSTOLM =   0.0
1577
        ELSEIF (ZARG1.GT.0.) THEN
1578
              LMSTOLAM =  90.0
1579
            ELSE
1580
              LMSTOLAM = -90.0
1581
            ENDIF
1582
      ELSE
1583
        LMSTOLM = ZRPI18*ATAN2(ZARG1,ZARG2)
1584
      ENDIF
1585
 
1586
      RETURN
1587
      END
1588
 
1589
 
1590
      REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
1591
C
1592
C**** PHSTOPH  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
1593
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
1594
C****                 ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1595
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1596
C**   AUFRUF   :   PHI = PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
1597
C**   ENTRIES  :   KEINE
1598
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
1599
C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
1600
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1601
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1602
C**   VERSIONS-
1603
C**   DATUM    :   03.05.90
1604
C**
1605
C**   EXTERNALS:   KEINE
1606
C**   EINGABE-
1607
C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
1608
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
1609
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
1610
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
1611
C**   AUSGABE-
1612
C**   PARAMETER:   WAHRE GEOGRAPHISCHE BREITE ALS WERT DER FUNKTION
1613
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
1614
C**
1615
C**   COMMON-
1616
C**   BLOECKE  :   KEINE
1617
C**
1618
C**   FEHLERBE-
1619
C**   HANDLUNG :   KEINE
1620
C**   VERFASSER:   D.MAJEWSKI
1621
 
1622
      REAL        LAMS,PHIS,POLPHI,POLLAM
1623
 
1624
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
1625
 
1626
      SINPOL = SIN(ZPIR18*POLPHI)
1627
      COSPOL = COS(ZPIR18*POLPHI)
1628
      ZPHIS  = ZPIR18*PHIS
1629
      ZLAMS  = LAMS
1630
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
1631
      ZLAMS  = ZPIR18*ZLAMS
1632
      ARG     = COSPOL*COS(ZPHIS)*COS(ZLAMS) + SINPOL*SIN(ZPHIS)
1633
 
1634
      PHSTOPH = ZRPI18*ASIN(ARG)
1635
 
1636
      RETURN
1637
      END
1638
 
1639
 
1640
      REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)
1641
C
1642
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
1643
C
1644
C**** LMTOLMS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM
1645
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
1646
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
1647
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1648
C**   AUFRUF   :   LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)
1649
C**   ENTRIES  :   KEINE
1650
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF
1651
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
1652
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1653
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1654
C**   VERSIONS-
1655
C**   DATUM    :   03.05.90
1656
C**
1657
C**   EXTERNALS:   KEINE
1658
C**   EINGABE-
1659
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
1660
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
1661
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
1662
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
1663
C**   AUSGABE-
1664
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
1665
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
1666
C**
1667
C**   COMMON-
1668
C**   BLOECKE  :   KEINE
1669
C**
1670
C**   FEHLERBE-
1671
C**   HANDLUNG :   KEINE
1672
C**   VERFASSER:   G. DE MORSIER
1673
 
1674
      REAL        LAM,PHI,POLPHI,POLLAM
1675
 
1676
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
1677
 
1678
      ZSINPOL = SIN(ZPIR18*POLPHI)
1679
      ZCOSPOL = COS(ZPIR18*POLPHI)
1680
      ZLAMPOL =     ZPIR18*POLLAM
1681
      ZPHI    =     ZPIR18*PHI
1682
      ZLAM    = LAM
1683
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
1684
      ZLAM    = ZPIR18*ZLAM
1685
 
1686
      ZARG1   = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI)
1687
      ZARG2   = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI)
1688
      IF (ABS(ZARG2).LT.1.E-30) THEN
1689
        IF (ABS(ZARG1).LT.1.E-30) THEN
1690
          LMTOLMS =   0.0
1691
        ELSEIF (ZARG1.GT.0.) THEN
1692
              LMTOLMS =  90.0
1693
            ELSE
1694
              LMTOLMS = -90.0
1695
            ENDIF
1696
      ELSE
1697
        LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
1698
      ENDIF
1699
 
1700
      RETURN
1701
      END
1702
 
1703
 
1704
      REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
1705
C
1706
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
1707
C
1708
C**** PHTOPHS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
1709
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
1710
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
1711
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1712
C**   AUFRUF   :   PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)
1713
C**   ENTRIES  :   KEINE
1714
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF
1715
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
1716
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1717
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1718
C**   VERSIONS-
1719
C**   DATUM    :   03.05.90
1720
C**
1721
C**   EXTERNALS:   KEINE
1722
C**   EINGABE-
1723
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
1724
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
1725
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
1726
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
1727
C**   AUSGABE-
1728
C**   PARAMETER:   ROTIERTE BREITE PHIS ALS WERT DER FUNKTION
1729
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
1730
C**
1731
C**   COMMON-
1732
C**   BLOECKE  :   KEINE
1733
C**
1734
C**   FEHLERBE-
1735
C**   HANDLUNG :   KEINE
1736
C**   VERFASSER:   G. DE MORSIER
1737
 
1738
      REAL        LAM,PHI,POLPHI,POLLAM
1739
 
1740
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
1741
 
1742
      ZSINPOL = SIN(ZPIR18*POLPHI)
1743
      ZCOSPOL = COS(ZPIR18*POLPHI)
1744
      ZLAMPOL = ZPIR18*POLLAM
1745
      ZPHI    = ZPIR18*PHI
1746
      ZLAM    = LAM
1747
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
1748
      ZLAM    = ZPIR18*ZLAM
1749
      ZARG    = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI)
1750
 
1751
      PHTOPHS = ZRPI18*ASIN(ZARG)
1752
 
1753
      RETURN
1754
      END