Subversion Repositories lagranto.ecmwf

Rev

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

Rev Author Line No. Line
3 michaesp 1
      PROGRAM trace
2
 
3
C     ********************************************************************
4
C     *                                                                  *
5
C     * Pseudo-lidar plots along trajectories                             *
6
C     *                                                                  *
7
C     * Heini Wernli       first version:       April 1993               *
8
C     * Michael Sprenger   major upgrade:       2008-2009                *
9
C     *                                                                  *
10
C     ********************************************************************
11
 
12
      implicit none
13
 
14
c     --------------------------------------------------------------------
15
c     Declaration of parameters
16
c     --------------------------------------------------------------------
17
 
18
c     Maximum number of levels for input files
19
      integer   nlevmax
20
      parameter (nlevmax=100)
21
 
22
c     Maximum number of input files (dates, length of trajectories)
23
      integer   ndatmax
24
      parameter (ndatmax=500)
25
 
26
c     Numerical epsilon (for float comparison)
27
      real      eps
28
      parameter (eps=0.001)
29
 
30
c     Conversion factors
31
      real      pi180                                   ! deg -> rad
32
      parameter (pi180=3.14159/180.)
33
      real      deg2km                                  ! deg -> km (at equator)
34
      parameter (deg2km=111.)
35
 
36
c     Prefix for primary and secondary fields
37
      character charp
38
      character chars
39
      parameter (charp='P')
40
      parameter (chars='S')
41
 
42
c     --------------------------------------------------------------------
43
c     Declaration of variables
44
c     --------------------------------------------------------------------
45
 
46
c     Input and output format for trajectories (see iotra.f)
47
      integer   inpmode
48
 
49
c     Input parameters
50
      character*80                           inpfile         ! Input trajectory file
51
      character*80                           outfile         ! Output netCDF file
52
      character*80                           outmode         ! Output mode (sum,mean)
53
      integer                                ntra            ! Number of trajectories
54
      integer                                ncol            ! Number of columns (including time, lon, lat, p)
55
      integer                                ntim            ! Number of times per trajectory
56
      integer                                ntrace0         ! Number of trace variables
57
      character*80                           tvar(200)       ! Tracing variable name (only the variable)
58
      character*1                            tfil(200)       ! Filename prefix 
59
      real                                   fac(200)        ! Scaling factor 
60
      integer                                compfl(200)     ! Computation flag (1=compute)
61
      integer                                numdat          ! Number of input files
62
      character*11                           dat(ndatmax)    ! Dates of input files
63
      real                                   timeinc         ! Time increment between input files
64
      real                                   tst             ! Time shift of start relative to first data file
65
      real                                   ten             ! Time shift of end relatiev to first data file  
66
      character*20                           startdate       ! First time/date on trajectory
67
      character*20                           enddate         ! Last time/date on trajectory
68
      character*80                           timecheck       ! Either 'yes' or 'no'
69
      character*80                           intmode         ! Interpolation mode ('normal', 'nearest')
70
	  real                                   pmin,pmax       ! Pressure range for output grid
71
	  integer                                npre            ! Number of pressure levels in output grid
72
	  character*80                           centering       ! Centering around trajectory position ('yes','no')
5 michaesp 73
	  character*80                       direction       ! Direction of lidar (vertical,lat,lon,normal)
74
      character*80                           dumpcoord       ! Dumping coordinates ('yes','no')
3 michaesp 75
 
76
c     Trajectories
5 michaesp 77
      real,allocatable, dimension (:,:,:) :: trainp                ! Input trajectories (ntra,ntim,ncol)
78
      integer                                reftime(6)            ! Reference date
79
      character*80                           varsinp(100)          ! Field names for input trajectory
80
      integer                                fid,fod               ! File identifier for inp and out trajectories
81
      real                                   x0_tra,y0_tra,p0_tra  ! Position of air parcel (physical space)
82
      real                                   reltpos0              ! Relative time of air parcel
83
      real                                   xind,yind,pind        ! Position of air parcel (grid space)
84
      integer                                fbflag                ! Flag for forward (1) or backward (-1) trajectories
3 michaesp 85
 
86
c     Meteorological fields from input file
87
      real,allocatable, dimension (:)     :: spt0,spt1       ! Surface pressure
88
      real,allocatable, dimension (:)     :: p3t0,p3t1       ! 3d-pressure 
89
      real,allocatable, dimension (:)     :: f3t0,f3t1       ! 3d field for tracing 
90
      character*80                           svars(100)      ! List of variables on S file
91
      character*80                           pvars(100)      ! List of variables on P file
92
      integer                                n_svars         ! Number of variables on S file
93
      integer                                n_pvars         ! Number of variables on P file
94
 
95
c     Input grid description
96
      real                                   pollon,pollat   ! Longitude/latitude of pole
97
      real                                   ak(100)         ! Vertical layers and levels
98
      real                                   bk(100) 
99
      real                                   xmin,xmax       ! Zonal grid extension
100
      real                                   ymin,ymax       ! Meridional grid extension
101
      integer                                nx,ny,nz        ! Grid dimensions
102
      real                                   dx,dy           ! Horizontal grid resolution
103
      integer                                hem             ! Flag for hemispheric domain
104
      integer                                per             ! Flag for periodic domain
105
      real                                   stagz           ! Vertical staggering
106
      real                                   mdv             ! Missing data value
107
 
108
c	  Output grid  and fields
109
      real                                   levels(1000)    ! Ouput levels
110
      real                                   times (1000)    ! Output times
111
	  real,allocatable, dimension (:,:)   :: out_pos         ! Position of trajectories
112
	  real,allocatable, dimension (:,:)   :: out_val         ! Output lidar field
113
	  real,allocatable, dimension (:,:)   :: out_cnt         ! # output lidar sum ups
114
 
115
 
116
c     Auxiliary variables
117
      integer                                i,j,k,l,n
118
      real                                   rd
119
      character*80                           filename
120
      real                                   time0,time1,reltpos
121
      integer                                itime0,itime1
122
      integer                                stat
123
      real                                   tstart
124
      integer                                iloaded0,iloaded1
125
      real                                   f0
126
      real                                   frac
127
      real                                   tload,tfrac
128
      integer                                isok
129
      character                              ch
130
      integer                                ind
131
      integer                                ind1,ind2,ind3,ind4,ind5
132
      integer                                ind6,ind7,ind8,ind9,ind0
133
      integer                                noutside
134
      real                                   delta
135
      integer                                itrace0
136
      character*80                           string
137
      character*80                           cdfname
138
      character*80                           varname
139
      real                                   time
140
      character*80                           longname
141
      character*80                           unit
142
      integer                                ind_time
143
      integer                                ind_pre
5 michaesp 144
      real                                   rlat,rlon
145
      real                                   x0,y0,p0
146
      real                                   vx0,vy0,vx1,vy1
147
      real                                   rotation,lon,lat
3 michaesp 148
 
149
c     Externals 
150
      real                                   int_index4
151
      external                               int_index4
152
 
153
c     --------------------------------------------------------------------
154
c     Start of program, Read parameters, get grid parameters
155
c     --------------------------------------------------------------------
156
 
157
c     Write start message
158
      print*,'========================================================='
159
      print*,'              *** START OF PROGRAM LIDAR ***'
160
      print*
161
 
162
c     Read parameters
163
      open(10,file='trace.param')
164
       read(10,*) inpfile
165
       read(10,*) outfile
166
       read(10,*) outmode
167
       read(10,*) startdate
168
       read(10,*) enddate 
169
       read(10,*) fbflag
170
       read(10,*) numdat
171
       if ( fbflag.eq.1) then
172
          do i=1,numdat
173
             read(10,'(a11)') dat(i)
174
          enddo
175
       else
176
          do i=numdat,1,-1
177
             read(10,'(a11)') dat(i)
178
          enddo
179
       endif
180
       read(10,*) timeinc
181
       read(10,*) tst
182
       read(10,*) ten
183
       read(10,*) ntra
184
       read(10,*) ntim
185
       read(10,*) ncol
186
       read(10,*) ntrace0
187
       do i=1,ntrace0
188
          read(10,*) tvar(i), fac(i), compfl(i), tfil(i)
189
       enddo
190
       read(10,*) n_pvars
191
       do i=1,n_pvars
192
          read(10,*) pvars(i)
193
       enddo
194
       read(10,*) n_svars
195
       do i=1,n_svars
196
          read(10,*) svars(i)
197
       enddo
198
       read(10,*) timecheck
199
       read(10,*) intmode
200
       read(10,*) pmin,pmax,npre
201
       read(10,*) centering
5 michaesp 202
       read(10,*) direction
203
       read(10,*) dumpcoord
3 michaesp 204
      close(10)
205
 
5 michaesp 206
c     Check that the direction is ok
207
      if ( ( direction.ne.'vertical' ).and.
208
     >     ( direction.ne.'lat'      ).and.
209
     >     ( direction.ne.'lon'      ).and.
210
     >     ( direction.ne.'normal'   ) ) 
211
     >then 
212
         print*,' ERROR: invalid direction ',trim(direction)
213
         stop
214
      endif
215
 
3 michaesp 216
c     Remove commented tracing fields
217
      itrace0 = 1
218
      do while ( itrace0.le.ntrace0) 
219
         string = tvar(itrace0)
220
         if ( string(1:1).eq.'#' ) then
221
            do i=itrace0,ntrace0-1
222
               tvar(i)   = tvar(i+1)
223
               fac(i)    = fac(i+1)
224
               compfl(i) = compfl(i+1)
225
               tfil(i)   = tfil(i+1)
226
            enddo
227
            ntrace0 = ntrace0 - 1
228
         else
229
            itrace0 = itrace0 + 1
230
         endif
231
      enddo
232
 
233
c     Set the formats of the input  files
234
      call mode_tra(inpmode,inpfile)
235
      if (inpmode.eq.-1) inpmode=1
236
 
237
C     Convert time shifts <tst,ten> from <hh.mm> into fractional time
238
      call hhmm2frac(tst,frac)
239
      tst = frac
240
      call hhmm2frac(ten,frac)
241
      ten = frac
242
 
243
c     Set the time for the first data file (depending on forward/backward mode)
244
      if (fbflag.eq.1) then
245
        tstart = -tst
246
      else
247
        tstart = tst
248
      endif
249
 
250
c     Read the constant grid parameters (nx,ny,nz,xmin,xmax,ymin,ymax,pollon,pollat)
251
c     The negative <-fid> of the file identifier is used as a flag for parameter retrieval  
252
      filename = charp//dat(1)
253
      varname  = 'U'
254
      call input_open (fid,filename)
255
      call input_grid (-fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
256
     >                 tstart,pollon,pollat,rd,rd,nz,rd,rd,rd,timecheck)
257
      call input_close(fid)
258
 
259
C     Allocate memory for some meteorological arrays
260
      allocate(spt0(nx*ny),stat=stat)
261
      if (stat.ne.0) print*,'*** error allocating array spt0 ***'   ! Surface pressure
262
      allocate(spt1(nx*ny),stat=stat)
263
      if (stat.ne.0) print*,'*** error allocating array spt1 ***'
264
      allocate(p3t0(nx*ny*nz),stat=stat)
265
      if (stat.ne.0) print*,'*** error allocating array p3t0 ***'   ! Pressure
266
      allocate(p3t1(nx*ny*nz),stat=stat)
267
      if (stat.ne.0) print*,'*** error allocating array p3t1 ***'
268
      allocate(f3t0(nx*ny*nz),stat=stat)
269
      if (stat.ne.0) print*,'*** error allocating array p3t0 ***'   ! Lidar field
270
      allocate(f3t1(nx*ny*nz),stat=stat)
271
      if (stat.ne.0) print*,'*** error allocating array p3t1 ***'
272
 
273
c	  Allocate memory for output field
274
	  allocate(out_pos(ntim,npre),stat=stat)
275
      if (stat.ne.0) print*,'*** error allocating array out_pos ***'
276
	  allocate(out_val(ntim,npre),stat=stat)
277
      if (stat.ne.0) print*,'*** error allocating array out_val ***'
278
	  allocate(out_cnt(ntim,npre),stat=stat)
279
      if (stat.ne.0) print*,'*** error allocating array out_cnt ***'
280
 
281
C     Get memory for trajectory arrays
282
      allocate(trainp(ntra,ntim,ncol),stat=stat)
283
      if (stat.ne.0) print*,'*** error allocating array tra      ***' 
284
 
285
c     Set the flags for periodic domains
286
      if ( abs(xmax-xmin-360.).lt.eps ) then
287
         per = 1
288
      elseif ( abs(xmax-xmin-360.+dx).lt.eps ) then
289
         per = 2
290
      else
291
         per = 0
292
      endif
293
 
294
C     Set logical flag for periodic data set (hemispheric or not)
295
      hem = 0
296
      if (per.eq.0.) then
297
         delta=xmax-xmin-360.
298
         if (abs(delta+dx).lt.eps) then               ! Program aborts: arrays must be closed
299
            print*,' ERROR: arrays must be closed... Stop'
300
         else if (abs(delta).lt.eps) then ! Periodic and hemispheric
301
           hem=1
302
           per=360.
303
        endif
304
      else                                            ! Periodic and hemispheric
305
         hem=1
306
      endif
307
 
308
c     Write some status information
309
      print*,'---- INPUT PARAMETERS -----------------------------------'
310
      print*
311
      print*,'  Input trajectory file  : ',trim(inpfile)
312
      print*,'  Format of input file   : ',inpmode
313
      print*,'  Output netCDF    file  : ',trim(outfile)
314
      print*,'  Format of output file  : ',trim(outmode)
315
      print*,'  Forward/backward       : ',fbflag
316
      print*,'  #tra                   : ',ntra
317
      print*,'  #col                   : ',ncol
318
      print*,'  #tim                   : ',ntim
319
      print*,'  No time check          : ',trim(timecheck)
320
      print*,'  Interpolation mode     : ',trim(intmode)
321
      do i=1,ntrace0
322
         if (compfl(i).eq.0) then
323
            print*,'  Tracing field          : ',
324
     >                 trim(tvar(i)), fac(i), ' 0 ', tfil(i)
325
         else
326
            print*,'  Tracing field          : ',
327
     >                trim(tvar(i)),' : online calc not supported'
328
         endif
329
      enddo
330
      print*,'  Output (pmin,pmax,n)   : ',pmin,pmax,npre
331
      print*,'  Centering              : ',trim(centering)
5 michaesp 332
      print*,'  Orientation            : ',trim(direction)
333
      print*,'  Coordinate Dump        : ',trim(dumpcoord)
3 michaesp 334
      print*
335
      print*,'---- INPUT DATA FILES -----------------------------------'
336
      print*
337
      call frac2hhmm(tstart,tload)
338
      print*,'  Time of 1st data file  : ',tload
339
      print*,'  #input files           : ',numdat
340
      print*,'  time increment         : ',timeinc
341
      call frac2hhmm(tst,tload)
342
      print*,'  Shift of start         : ',tload
343
      call frac2hhmm(ten,tload)
344
      print*,'  Shift of end           : ',tload
345
      print*,'  First/last input file  : ',trim(dat(1)),
346
     >                                     ' ... ',
5 michaesp 347
     >                                     trim(dat(numdat)) 
3 michaesp 348
      print*,'  Primary variables      : ',trim(pvars(1))
349
      do i=2,n_pvars
350
         print*,'                         : ',trim(pvars(i))
351
      enddo
352
      if ( n_svars.ge.1 ) then
353
         print*,'  Secondary variables    : ',trim(svars(1))
354
         do i=2,n_svars
355
            print*,'                         : ',trim(svars(i))
356
         enddo
357
      endif
358
      print*
359
      print*,'---- CONSTANT GRID PARAMETERS ---------------------------'
360
      print*
361
      print*,'  xmin,xmax     : ',xmin,xmax
362
      print*,'  ymin,ymax     : ',ymin,ymax
363
      print*,'  dx,dy         : ',dx,dy
364
      print*,'  pollon,pollat : ',pollon,pollat
365
      print*,'  nx,ny,nz      : ',nx,ny,nz
366
      print*,'  per, hem      : ',per,hem
367
      print*
368
 
369
c     --------------------------------------------------------------------
370
c     Load the input trajectories
371
c     --------------------------------------------------------------------
372
 
373
c     Read the input trajectory file
374
      call ropen_tra(fid,inpfile,ntra,ntim,ncol,reftime,varsinp,inpmode)
375
      call read_tra (fid,trainp,ntra,ntim,ncol,inpmode)
376
      call close_tra(fid,inpmode)
377
 
378
c     Check that first four columns correspond to time,lon,lat,p
379
      if ( (varsinp(1).ne.'time' ).or.
380
     >     (varsinp(2).ne.'xpos' ).and.(varsinp(2).ne.'lon' ).or.
381
     >     (varsinp(3).ne.'ypos' ).and.(varsinp(3).ne.'lat' ).or.
382
     >     (varsinp(4).ne.'ppos' ).and.(varsinp(4).ne.'p'   ) )
383
     >then
384
         print*,' ERROR: problem with input trajectories ...'
385
         stop
386
      endif
387
      varsinp(1) = 'time'
388
      varsinp(2) = 'lon'
389
      varsinp(3) = 'lat'
390
      varsinp(4) = 'p'
391
 
392
c     Write some status information of the input trajectories
393
      print*,'---- INPUT TRAJECTORIES ---------------------------------'
394
      print*
395
      print*,' Start date             : ',trim(startdate)
396
      print*,' End date               : ',trim(enddate)
397
      print*,' Reference time (year)  : ',reftime(1)
398
      print*,'                (month) : ',reftime(2)
399
      print*,'                (day)   : ',reftime(3)
400
      print*,'                (hour)  : ',reftime(4)
401
      print*,'                (min)   : ',reftime(5)
402
      print*,' Time range (min)       : ',reftime(6)
403
      do i=1,ncol
404
         print*,' Var                    :',i,trim(varsinp(i))
405
      enddo
406
      print*
407
 
408
c     Check that first time is 0 - otherwise the tracing will produce
409
c     wrong results because later in the code only absolute times are
410
c     considered: <itime0   = int(abs(tfrac-tstart)/timeinc) + 1>. This 
411
c     will be changed in a future version.
412
      if ( abs( trainp(1,1,1) ).gt.eps ) then
413
         print*,' ERROR: First time of trajectory must be 0, i.e. '
414
         print*,'     correspond to the reference date. Otherwise'
415
         print*,'     the tracing will give wrong results... STOP'
416
         stop
417
      endif
5 michaesp 418
 
419
c     If requested, open the coordinate dump file
420
      if ( dumpcoord.eq.'yes' ) then
421
         open(10,file=trim(outfile)//'.coord')
422
      endif
423
 
3 michaesp 424
c     --------------------------------------------------------------------
425
c     Trace the fields (fields available on input files)
426
c     --------------------------------------------------------------------
427
 
428
      print*
429
      print*,'---- LIDAR FROM PRIMARY AND SECONDARY DATA FILES ------'
430
 
431
c     Loop over all tracing fields
432
      do i=1,ntrace0
433
 
434
c	      Skip all fields marked for online calculation
435
          if ( compfl(i).eq.1 ) goto 110
436
 
437
c	      Init the output fields: position and lidar field
438
	      do k=1,ntim
439
	      	do l=1,npre
440
	      		out_pos(k,l) = 0.
441
	      	    out_val(k,l) = 0.
442
	      	    out_cnt(k,l) = 0.
443
	      	 enddo
444
	      enddo
445
 
446
c         Write some status information
447
          print*
448
          print*,' Now lidaring           : ',
449
     >         trim(tvar(i)),compfl(i),' ',trim(tfil(i))
450
 
451
c         Reset flags for load manager
452
          iloaded0 = -1
453
          iloaded1 = -1
454
 
455
c         Reset the counter for fields outside domain
456
          noutside = 0
457
 
458
c         Loop over all times
459
          do j=1,ntim
460
 
461
c            Convert trajectory time from hh.mm to fractional time
462
             call hhmm2frac(trainp(1,j,1),tfrac)
463
 
464
c            Get the times which are needed
465
             itime0   = int(abs(tfrac-tstart)/timeinc) + 1
466
             time0    = tstart + fbflag * real(itime0-1) * timeinc
467
             itime1   = itime0 + 1
468
             time1    = time0 + fbflag * timeinc
469
             if ( itime1.gt.numdat ) then
470
                itime1 = itime0
471
                time1  = time0
472
             endif
473
 
474
c            Load manager: Check whether itime0 can be copied from itime1
475
             if ( itime0.eq.iloaded1 ) then
476
                f3t0     = f3t1
477
                p3t0     = p3t1
478
                spt0     = spt1
479
                iloaded0 = itime0
480
             endif
481
 
482
c            Load manager: Check whether itime1 can be copied from itime0
483
             if ( itime1.eq.iloaded0 ) then
484
                f3t1     = f3t0
485
                p3t1     = p3t0
486
                spt1     = spt0
487
                iloaded1 = itime1
488
             endif
489
 
490
c            Load manager:  Load first time (tracing variable and grid)
491
             if ( itime0.ne.iloaded0 ) then
492
 
493
                filename = tfil(i)//dat(itime0)
494
                call frac2hhmm(time0,tload)
495
                varname  = tvar(i) 
496
                write(*,'(a23,a20,a3,a5,f7.2)') 
497
     >               '    ->  loading          : ',
498
     >               trim(filename),'  ',trim(varname),tload
499
                call input_open (fid,filename)                
500
                call input_wind 
501
     >               (fid,varname,f3t0,tload,stagz,mdv,
502
     >               xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck) 
503
 
504
                call input_grid      
505
     >               (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
506
     >                tload,pollon,pollat,p3t0,spt0,nz,ak,bk,stagz,
507
     >                timecheck)
508
                call input_close(fid)
509
 
510
                iloaded0 = itime0
511
 
512
             endif
513
 
514
c            Load manager: Load second time (tracing variable and grid)
515
             if ( itime1.ne.iloaded1 ) then
516
 
517
                filename = tfil(i)//dat(itime1)
518
                call frac2hhmm(time1,tload)
519
                varname  = tvar(i) 
520
                write(*,'(a23,a20,a3,a5,f7.2)') 
521
     >               '    ->  loading          : ',
522
     >               trim(filename),'  ',trim(varname),tload
523
                call input_open (fid,filename)
524
                call input_wind 
525
     >               (fid,varname,f3t1,tload,stagz,mdv,
526
     >               xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)      
527
                call input_grid      
528
     >               (fid,varname,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
529
     >               tload,pollon,pollat,p3t1,spt1,nz,ak,bk,stagz,
530
     >               timecheck)
531
                call input_close(fid)
532
 
533
                iloaded1 = itime1
534
 
535
             endif
536
 
537
c            Loop over all trajectories
538
             do k=1,ntra
539
 
5 michaesp 540
c               Set the trajectory position 
541
                x0_tra = trainp(k,j,2)                          ! Longitude
542
                y0_tra = trainp(k,j,3)                          ! Latitude
543
                p0_tra = trainp(k,j,4)                          ! Pressure
3 michaesp 544
 
5 michaesp 545
c               Get rotation angle - orient normal to trajectory
546
                if ( direction.eq.'normal' ) then
547
 
548
                   vx0 = 1.                     
549
                   vy0 = 0.
550
 
551
                   if ( j.lt.ntim ) then
552
                      lat =  0.5 * ( trainp(k,j,3) + trainp(k,j+1,3) )
553
                      vx1 = ( trainp(k,j+1,2) - trainp(k,j,2) ) * 
554
     >                      cos( lat * pi180 )
555
                      vy1 = ( trainp(k,j+1,3) - trainp(k,j,3) )
556
                   else
557
                      lat =  0.5 * ( trainp(k,j,3) + trainp(k,j-1,3) )
558
                      vx1 = ( trainp(k,j,2) - trainp(k,j-1,2) ) * 
559
     >                      cos( lat * pi180 )
560
                      vy1 = ( trainp(k,j,3) - trainp(k,j-1,3) )
561
                   endif
562
 
563
                   if ( vx1.gt.180  ) vx1 = vx1 - 360
564
                   if ( vx1.lt.-180 ) vx1 = vx1 + 360.
565
 
566
                   call getangle (vx0,vy0,vx1,vy1,rotation)
567
                   rotation = -rotation
568
 
569
                else
570
                   rotation = 0.
571
                endif
572
 
3 michaesp 573
c               Set the relative time
574
                call hhmm2frac(trainp(k,j,1),tfrac)
575
                reltpos0 = fbflag * (tfrac-time0)/timeinc    
576
 
5 michaesp 577
c               Loop over pressure profile (or other positions for horizontal mode)
578
	        do l=1,npre
3 michaesp 579
 
5 michaesp 580
c                     Vertical
581
                      if ( direction.eq.'vertical' ) then
582
                        x0 = x0_tra
583
                        y0 = y0_tra
584
                        p0 = pmin + real(l-1)/real(npre-1) * (pmax-pmin)
585
                        if ( centering.eq.'yes' )then
586
                           p0 = p0 + trainp(k,j,4)
587
                        endif
3 michaesp 588
 
5 michaesp 589
c                     Longitude
590
                      elseif ( direction.eq.'lon' ) then
591
                        x0 = pmin + real(l-1)/real(npre-1) * (pmax-pmin)
592
                        y0 = y0_tra
593
                        p0 = p0_tra
594
                        if ( centering.eq.'yes' )then
595
                           x0 = x0 + x0_tra
596
                        endif
597
 
598
c                     Latitude
599
                      elseif ( direction.eq.'lat' ) then
600
                        x0 = x0_tra
601
                        y0 = pmin + real(l-1)/real(npre-1) * (pmax-pmin)
602
                        p0 = p0_tra
603
                        if ( centering.eq.'yes' )then
604
                           y0 = y0 + y0_tra
605
                        endif
3 michaesp 606
 
5 michaesp 607
c                     Normal to trajerctory
608
                      elseif ( direction.eq.'normal' ) then
609
 
610
c                        Set the coordinate in the rotated system
611
                         rlat = pmin + 
612
     >                             real(l-1)/real(npre-1) * (pmax-pmin)
613
                         rlon = 0.
614
 
615
c                        Transform it back to geographical lon/lat
616
                         call getenvir_b (x0_tra,y0_tra,rotation,
617
     >                                              x0,y0,rlon,rlat,1)
618
 
619
c                        Pressure unchanged
620
                         p0 = p0_tra
621
 
622
                      endif
623
 
624
c                     Handle periodic boundaries in zonal direction
625
                      if ( (x0.gt.xmax).and.(per.ne.0) ) x0 = x0 - 360.
626
                      if ( (x0.lt.xmin).and.(per.ne.0) ) x0 = x0 + 360.
627
 
628
c                     Handle pole problems for hemispheric data (taken from caltra.f)
629
                      if ((hem.eq.1).and.(y0.gt.90.)) then
630
                         y0=180.-y0
631
                         x0=x0+per/2.
632
                      endif
633
                      if ((hem.eq.1).and.(y0.lt.-90.)) then
634
                         y0=-180.-y0
635
                         x0=x0+per/2.
636
                      endif
637
                      if (y0.gt.89.99) then
638
                         y0=89.99
639
                      endif     
640
 
641
c                 If requested, dump the lidar coordinates
642
                  if ( (dumpcoord.eq.'yes').and.(i.eq.1) ) then
643
                     write(10,'3f10.2') x0,y0,trainp(k,j,1)
644
                     write(10,'3f10.2') x0_tra,y0_tra,5.
645
                  endif
646
 
3 michaesp 647
C                 Get the index where to interpolate (x0,y0,p0)
648
                  if ( (abs(x0-mdv).gt.eps).and.
649
     >                 (abs(y0-mdv).gt.eps) )
650
     >            then
651
                     call get_index4 (xind,yind,pind,x0,y0,p0,reltpos0,
652
     >                                p3t0,p3t1,spt0,spt1,3,
653
     >                                nx,ny,nz,xmin,ymin,dx,dy,mdv)
654
                  else
655
                     xind = mdv
656
                     yind = mdv
657
                     pind = mdv
658
                  endif
659
 
660
c                 If requested, apply nearest-neighbor interpolation
661
                  if ( intmode.eq.'nearest') then
662
 
663
                     xind = real( nint(xind) )
664
                     yind = real( nint(yind) )
665
                     pind = real( nint(pind) )
666
 
667
                     if ( xind.lt.1.  ) xind = 1.
668
                     if ( xind.gt.nx  ) xind = real(nx)
669
                     if ( yind.lt.1.  ) yind = 1.
670
                     if ( yind.gt.ny  ) yind = real(ny)
671
 
672
                     if ( pind.lt.1.  ) pind = 1.
673
                     if ( pind.gt.nz  ) pind = real(nz)
674
 
675
                  endif
676
 
677
c                 Do the interpolation: everthing is ok
678
                  if ( (xind.ge.1.).and.(xind.le.real(nx)).and.
679
     >                 (yind.ge.1.).and.(yind.le.real(ny)).and.
680
     >                 (pind.ge.1.).and.(pind.le.real(nz)) )
681
     >            then
682
                     f0 = int_index4(f3t0,f3t1,nx,ny,nz,
683
     >                               xind,yind,pind,reltpos0,mdv)
684
 
685
c                 Set to missing data
686
                  else
687
                     f0       = mdv
688
                  endif
689
 
690
c	              Save result to output array
691
                  if (abs(f0-mdv).gt.eps) then
692
                     out_val(j,l) = out_val(j,l) + f0 * fac(i)
693
                     out_cnt(j,l) = out_cnt(j,l) + 1.
694
 
695
	              endif
696
 
697
c              End loop over all pressure levels
698
	           enddo
699
 
5 michaesp 700
c	           Save output - time index
3 michaesp 701
	           ind_time = j
5 michaesp 702
 
703
c                  Save output - space index for 'no centering'
3 michaesp 704
	           if ( centering.eq.'no' ) then
5 michaesp 705
                      if ( direction.eq.'vertical') then 
706
                         ind_pre  = nint( real(npre) *
707
     >          	       ( (p0_tra - pmin)/(pmax-pmin) ) + 1.)
708
                      elseif ( direction.eq.'lon') then 
709
                         ind_pre  = nint( real(npre) *
710
     >          	       ( (x0_tra - pmin)/(pmax-pmin) ) + 1.)
711
                      elseif ( direction.eq.'lat') then 
712
                         ind_pre  = nint( real(npre) *
713
     >          	       ( (y0_tra - pmin)/(pmax-pmin) ) + 1.)
714
                      endif
715
 
716
c                  Save output - space index for 'centering'
3 michaesp 717
	           else
718
	           	  ind_pre  = nint( real(npre) *
719
     >          	       ( (0.            - pmin)/(pmax-pmin) ) + 1.)
720
	           endif
721
 
5 michaesp 722
c                  Update the output array
3 michaesp 723
	           if ( (ind_time.ge.1).and.(ind_time.le.ntim).and.
724
     >              (ind_pre .ge.1).and.(ind_pre .le.npre) )
725
     >         then
726
                    out_pos(ind_time,ind_pre) =
727
     >                          	out_pos(ind_time,ind_pre) + 1.
728
	           endif
729
 
730
c	         End loop over all trajectories
731
             enddo
732
 
733
c	      End loop over all times
734
          enddo
735
 
736
c	      Write the trajectory position to netCDF file - only once
737
	      if ( i.eq.1 ) then
738
	      	  cdfname  = outfile
739
	      	  varname  = 'POSITION'
740
	      	  longname = 'position of trajectory points'
741
	      	  unit     = 'none'
742
	      	  time     = 0.
743
              do k=1,npre
744
              	levels(k) = pmin + real(k-1)/real(npre-1) * (pmax-pmin)
745
              enddo
746
              do k=1,ntim
747
                 times(k) = trainp(1,k,1)
748
              enddo
749
              call writecdf2D_cf
750
     >            (cdfname,varname,longname,unit,out_pos,time,levels,
5 michaesp 751
     >             times,npre,ntim,1,1,direction)
3 michaesp 752
	      endif
753
 
754
c	      If no valid lidar count: set the field to missing data
755
          do k=1,ntim
756
          	do l=1,npre
757
          		if (abs(out_cnt(k,l)).lt.eps) then
758
          			out_val(k,l) = mdv
759
          	    endif
760
          	 enddo
761
          enddo
762
 
763
c	      If requested, calculate the mean of the lidar field
764
	      if ( outmode.eq.'mean' ) then
765
	      	do k=1,ntim
766
	      		do l=1,npre
767
	      			if ( (abs(out_val(k,l)-mdv).gt.eps).and.
768
     >                   (abs(out_cnt(k,l)    ).gt.0. ) )
769
     >              then
770
	      				out_val(k,l) = out_val(k,l) / out_cnt(k,l)
771
	      		    endif
772
	      		 enddo
773
	          enddo
774
	      endif
775
 
776
c	      Write the lidar field and count
777
	      cdfname  = outfile
778
	      if (outmode.eq.'sum' ) then
779
	         varname  = trim(tvar(i))//'_SUM'
780
	      elseif (outmode.eq.'mean' ) then
781
	         varname  = trim(tvar(i))//'_MEAN'
782
	      endif
783
	      longname = 'sum over all '//trim(tvar(i))//' profiles'
784
	      unit     = 'not given'
785
	      time     = 0.
786
          call writecdf2D_cf
787
     >            (cdfname,varname,longname,unit,out_val,time,levels,
5 michaesp 788
     >             times,npre,ntim,0,1,direction)
3 michaesp 789
 
790
	  	  cdfname  = outfile
791
	      varname  = trim(tvar(i))//'_CNT'
792
	      longname = 'counts of all '//trim(tvar(i))//' profiles'
793
	      unit     = 'not given'
794
	      time     = 0.
795
          call writecdf2D_cf
796
     >            (cdfname,varname,longname,unit,out_cnt,time,levels,
5 michaesp 797
     >             times,npre,ntim,0,1,direction)
3 michaesp 798
 
799
c         Exit point for loop over all tracing variables
800
 110      continue
801
 
802
c	   End loop over all lidar variables
803
       enddo
804
 
805
 
806
c     --------------------------------------------------------------------
807
c     Write output to netCDF file
808
c     --------------------------------------------------------------------
809
 
810
c     Write status information
811
      print*
812
      print*,'---- WRITE OUTPUT LIDAR FIELDS --------------------------'
813
      print*
814
 
5 michaesp 815
c     Close coord dump file
816
      print*,' LIDAR written to      : ',trim(outfile)
817
      if ( dumpcoord.eq.'yes' ) then
818
        print*,' Coordinates dumped to : ',trim(outfile)//'.coord'
819
      endif
3 michaesp 820
 
821
c     Write some status information, and end of program message
822
      print*  
823
      print*,'---- STATUS INFORMATION --------------------------------'
824
      print*
825
      print*,' ok'
826
      print*
827
      print*,'              *** END OF PROGRAM LIDAR ***'
828
      print*,'========================================================='
829
 
830
 
831
      end 
832
 
833
 
834
c     ********************************************************************
835
c     * INPUT / OUTPUT SUBROUTINES                                       *
836
c     ********************************************************************
837
 
838
c     --------------------------------------------------------------------
839
c     Subroutines to write 2D CF netcdf output file
840
c     --------------------------------------------------------------------
841
 
842
      subroutine writecdf2D_cf
843
     >          (cdfname,varname,longname,unit,arr,time,levels,times,
5 michaesp 844
     >           npre,ntim,crefile,crevar,direction)
3 michaesp 845
 
846
c     Create and write to the CF netcdf file <cdfname>. The variable
847
c     with name <varname> and with time <time> is written. The data
848
c     are in the two-dimensional array <arr>.  The flags <crefile> and
849
c     <crevar> determine whether the file and/or the variable should
850
c     be created.
851
 
852
      USE netcdf
853
 
854
      IMPLICIT NONE
855
 
856
c     Declaration of input parameters
857
      character*80 cdfname
858
      character*80 varname,longname,unit
859
      integer      npre,ntim
860
      real         arr(ntim,npre)
861
      real         levels(npre)
862
      real         times (ntim)
863
      real         time
864
      integer      crefile,crevar
5 michaesp 865
      character*80 direction
3 michaesp 866
 
867
c     Numerical epsilon
868
      real         eps
869
      parameter    (eps=1.e-5)
870
 
871
c     Local variables
872
      integer      ierr
873
      integer      ncID
874
      integer      LevDimId,    varLevID
875
      integer      TimeDimID,   varTimeID
876
      real         timeindex
877
      integer      i,j
878
      integer      nvars,varids(100)
879
      integer      ndims,dimids(100)
880
      real         timelist(1000)
881
      integer      ntimes
882
      integer      ind
883
      integer      varID
884
 
885
c     Quick an dirty solution for fieldname conflict
886
      if ( varname.eq.'time' ) varname = 'TIME'
887
 
888
c     Initially set error to indicate no errors.
889
      ierr = 0
890
 
891
c     ---- Create the netCDF - skip if <crefile=0> ----------------------
892
      if ( crefile.ne.1 ) goto 100
893
 
894
c     Create the file
895
      ierr = nf90_create(trim(cdfname), NF90_CLOBBER, ncID)
896
 
897
c     Define dimensions
898
      ierr=nf90_def_dim(ncID,'level',npre, LevDimID )
899
      ierr=nf90_def_dim(ncID,'time' ,ntim, TimeDimID)
900
 
5 michaesp 901
c     Define space coordinate 
3 michaesp 902
      ierr = nf90_def_var(ncID,'level',NF90_FLOAT,
903
     >     (/ LevDimID /),varLevID)
5 michaesp 904
      if ( direction.eq.'vertical' ) then
905
        ierr = nf90_put_att(ncID, varLevID, "standard_name","level")
906
        ierr = nf90_put_att(ncID, varLevID, "units"        ,"hPa")
907
      elseif ( direction.eq.'lat' ) then
908
        ierr = nf90_put_att(ncID, varLevID, "standard_name","latitude")
909
        ierr = nf90_put_att(ncID, varLevID, "units"        ,"deg")
910
      elseif ( direction.eq.'lon' ) then
911
        ierr = nf90_put_att(ncID, varLevID, "standard_name","longitude")
912
        ierr = nf90_put_att(ncID, varLevID, "units"        ,"deg")
913
      elseif ( direction.eq.'normal' ) then
914
        ierr = nf90_put_att(ncID, varLevID, "standard_name","normal")
915
        ierr = nf90_put_att(ncID, varLevID, "units"        ,"deg")
916
      endif
3 michaesp 917
 
5 michaesp 918
c     Define time coordinate
3 michaesp 919
      ierr = nf90_def_var(ncID,'time',NF90_FLOAT,
920
     >     (/ TimeDimID /), varTimeID)
921
      ierr = nf90_put_att(ncID, varTimeID, "long_name",    "time")
922
      ierr = nf90_put_att(ncID, varTimeID, "units",       "hours")
923
 
924
c     Write global attributes
925
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'Conventions', 'CF-1.0')
926
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'title',
927
     >     'pseudo-lidar from trajectory file')
928
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'source',
929
     >     'Lagranto Trajectories')
930
      ierr = nf90_put_att(ncID, NF90_GLOBAL, 'institution',
931
     >     'ETH Zurich, IACETH')
932
 
933
c     Check whether the definition was successful
934
      ierr = nf90_enddef(ncID)
935
      if (ierr.gt.0) then
936
         print*, 'An error occurred while attempting to ',
937
     >        'finish definition mode.'
938
         stop
939
      endif
940
 
941
c     Write coordinate data
942
      ierr = nf90_put_var(ncID,varLevID  ,levels)
943
      ierr = nf90_put_var(ncID,varTimeID ,times )
944
 
945
c     Close netCDF file
946
      ierr = nf90_close(ncID)
947
 
948
 100  continue
949
 
950
c     ---- Define a new variable - skip if <crevar=0> -----------------------
951
 
952
      if ( crevar.ne.1 ) goto 110
953
 
5 michaesp 954
      print*,'Now defining ',trim(varname)
955
 
3 michaesp 956
c     Open the file for read/write access
957
      ierr = nf90_open  (trim(cdfname), NF90_WRITE  , ncID)
958
 
959
c     Get the IDs for dimensions
960
      ierr = nf90_inq_dimid(ncID,'level', LevDimID )
961
      ierr = nf90_inq_dimid(ncID,'time' , TimeDimID)
962
 
963
c     Enter define mode
964
      ierr = nf90_redef(ncID)
965
 
966
c     Write definition and add attributes
967
      ierr = nf90_def_var(ncID,varname,NF90_FLOAT,
5 michaesp 968
     >                   (/ TimeDimID, LevDimID /),varID)
3 michaesp 969
      ierr = nf90_put_att(ncID, varID, "long_name" , longname )
970
      ierr = nf90_put_att(ncID, varID, "units"     , unit     )
971
      ierr = nf90_put_att(ncID, varID, '_FillValue', -999.99  )
972
 
973
c     Check whether definition was successful
974
      ierr = nf90_enddef(ncID)
975
      if (ierr.gt.0) then
976
         print*, 'An error occurred while attempting to ',
977
     >           'finish definition mode.'
978
         stop
979
      endif
980
 
981
c     Close netCDF file
982
      ierr = nf90_close(ncID)
983
 
984
 110  continue
985
 
986
c     ---- Write data --------------------------------------------------
987
 
988
c     Open the file for read/write access
989
      ierr = nf90_open  (trim(cdfname), NF90_WRITE , ncID)
990
 
991
c     Get the varID
992
      ierr = nf90_inq_varid(ncID,varname, varID )
993
      if (ierr.ne.0) then
994
         print*,'Variable ',trim(varname),' is not defined on ',
995
     >          trim(cdfname)
996
         stop
997
      endif
998
 
999
c     Write data block
1000
      ierr = nf90_put_var(ncID,varID,arr,
1001
     >                    start = (/ 1, 1 /),
1002
     >                    count = (/ ntim, npre/) )
1003
 
1004
c     Check whether writing was successful
1005
      ierr = nf90_close(ncID)
1006
      if (ierr.ne.0) then
1007
         write(*,*) trim(nf90_strerror(ierr))
1008
         write(*,*) 'An error occurred while attempting to ',
1009
     >              'close the netcdf file.'
1010
         write(*,*) 'in clscdf_CF'
1011
      endif
1012
 
1013
      end
1014
 
5 michaesp 1015
c     ********************************************************************************
1016
c     * Coordinate rotation - lidar normal to trajectory                             *
1017
c     ********************************************************************************
3 michaesp 1018
 
5 michaesp 1019
c     --------------------------------------------------------------------------------
1020
c     Backward coordinate transformation (Rotated lon/lat -> True lon/lat)
1021
c     --------------------------------------------------------------------------------
3 michaesp 1022
 
5 michaesp 1023
      SUBROUTINE getenvir_b (clon,clat,rotation,
1024
     >                       lon,lat,rlon,rlat,n)
1025
 
1026
      implicit none
1027
 
1028
c     Declaration of input and output parameters
1029
      integer     n
1030
      real        clon,clat,rotation
1031
      real        lon(n), lat(n)
1032
      real        rlon(n),rlat(n)
1033
 
1034
c     Auxiliary variables 
1035
      real         pollon,pollat
1036
      integer      i
1037
      real         rlon1(n),rlat1(n)
1038
 
1039
c     Externals
1040
      real         lmstolm,phstoph
1041
      external     lmstolm,phstoph
1042
 
1043
c     First coordinate transformation (make the local coordinate system parallel to equator)
1044
      pollon=-180.
1045
      pollat=90.+rotation
1046
      do i=1,n
1047
         rlon1(i)=90.+lmstolm(rlat(i),rlon(i)-90.,pollat,pollon)
1048
         rlat1(i)=phstoph(rlat(i),rlon(i)-90.,pollat,pollon)            
1049
      enddo
1050
 
1051
c     Second coordinate transformation (make the local coordinate system parallel to equator)
1052
      pollon=clon-180.
1053
      if (pollon.lt.-180.) pollon=pollon+360.
1054
      pollat=90.-clat
1055
      do i=1,n
1056
         lon(i)=lmstolm(rlat1(i),rlon1(i),pollat,pollon)
1057
         lat(i)=phstoph(rlat1(i),rlon1(i),pollat,pollon)            
1058
      enddo
1059
 
1060
      END
1061
 
1062
c     ---------------------------------------------------------------------
1063
c     Determine the angle between two vectors
1064
c     ---------------------------------------------------------------------
1065
 
1066
      SUBROUTINE getangle (vx1,vy1,vx2,vy2,angle)
1067
 
1068
c     Given two vectors <vx1,vy1> and <vx2,vy2>, determine the angle (in deg) 
1069
c     between the two vectors.
1070
 
1071
      implicit none
1072
 
1073
c     Declaration of subroutine parameters
1074
      real vx1,vy1
1075
      real vx2,vy2
1076
      real angle
1077
 
1078
c     Auxiliary variables and parameters
1079
      real len1,len2,len3
1080
      real val1,val2,val3
1081
      real pi
1082
      parameter (pi=3.14159265359)
1083
 
1084
      len1=sqrt(vx1*vx1+vy1*vy1)
1085
      len2=sqrt(vx2*vx2+vy2*vy2)
1086
 
1087
      if ((len1.gt.0.).and.(len2.gt.0.)) then
1088
         vx1=vx1/len1
1089
         vy1=vy1/len1
1090
         vx2=vx2/len2
1091
         vy2=vy2/len2
3 michaesp 1092
 
5 michaesp 1093
         val1=vx1*vx2+vy1*vy2
1094
         val2=-vy1*vx2+vx1*vy2
3 michaesp 1095
 
5 michaesp 1096
         len3=sqrt(val1*val1+val2*val2)
1097
 
1098
         if ( (val1.ge.0.).and.(val2.ge.0.) ) then
1099
            val3=acos(val1/len3)
1100
         else if ( (val1.lt.0.).and.(val2.ge.0.) ) then
1101
            val3=pi-acos(abs(val1)/len3)
1102
         else if ( (val1.ge.0.).and.(val2.le.0.) ) then
1103
            val3=-acos(val1/len3)
1104
         else if ( (val1.lt.0.).and.(val2.le.0.) ) then
1105
            val3=-pi+acos(abs(val1)/len3)
1106
         endif
1107
      else
1108
         val3=0.
1109
      endif
3 michaesp 1110
 
5 michaesp 1111
      angle=180./pi*val3                                                                                     
3 michaesp 1112
 
5 michaesp 1113
      END
3 michaesp 1114
 
5 michaesp 1115
c     --------------------------------------------------------------------------------
1116
c     Transformation routine: LMSTOLM and PHSTOPH from library gm2em            
1117
c     --------------------------------------------------------------------------------
1118
 
1119
      REAL FUNCTION LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
1120
C
1121
C**** LMSTOLM  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
1122
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
1123
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
1124
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1125
C**   AUFRUF   :   LAM = LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
1126
C**   ENTRIES  :   KEINE
1127
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
1128
C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
1129
C**                IM ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1130
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1131
C**   VERSIONS-
1132
C**   DATUM    :   03.05.90
1133
C**
1134
C**   EXTERNALS:   KEINE
1135
C**   EINGABE-
1136
C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
1137
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
1138
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
1139
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
1140
C**   AUSGABE-
1141
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
1142
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
1143
C**
1144
C**   COMMON-
1145
C**   BLOECKE  :   KEINE
1146
C**
1147
C**   FEHLERBE-
1148
C**   HANDLUNG :   KEINE
1149
C**   VERFASSER:   D.MAJEWSKI
1150
 
1151
      REAL        LAMS,PHIS,POLPHI,POLLAM
1152
 
1153
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
1154
 
1155
      ZSINPOL = SIN(ZPIR18*POLPHI)
1156
      ZCOSPOL = COS(ZPIR18*POLPHI)
1157
      ZLAMPOL = ZPIR18*POLLAM
1158
      ZPHIS   = ZPIR18*PHIS
1159
      ZLAMS   = LAMS
1160
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
1161
      ZLAMS   = ZPIR18*ZLAMS
1162
 
1163
      ZARG1   = SIN(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
1164
     1                          ZCOSPOL*           SIN(ZPHIS)) -
1165
     2          COS(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
1166
      ZARG2   = COS(ZLAMPOL)*(- ZSINPOL*COS(ZLAMS)*COS(ZPHIS)  +
1167
     1                          ZCOSPOL*           SIN(ZPHIS)) +
1168
     2          SIN(ZLAMPOL)*           SIN(ZLAMS)*COS(ZPHIS)
1169
      IF (ABS(ZARG2).LT.1.E-30) THEN
1170
        IF (ABS(ZARG1).LT.1.E-30) THEN
1171
          LMSTOLM =   0.0
1172
        ELSEIF (ZARG1.GT.0.) THEN
1173
              LMSTOLAM =  90.0
1174
            ELSE
1175
              LMSTOLAM = -90.0
1176
            ENDIF
1177
      ELSE
1178
        LMSTOLM = ZRPI18*ATAN2(ZARG1,ZARG2)
1179
      ENDIF
1180
 
1181
      RETURN
1182
      END
1183
 
1184
 
1185
      REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
1186
C
1187
C**** PHSTOPH  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
1188
C****                 EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
1189
C****                 ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1190
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1191
C**   AUFRUF   :   PHI = PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
1192
C**   ENTRIES  :   KEINE
1193
C**   ZWECK    :   BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
1194
C**                EINEN PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
1195
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1196
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1197
C**   VERSIONS-
1198
C**   DATUM    :   03.05.90
1199
C**
1200
C**   EXTERNALS:   KEINE
1201
C**   EINGABE-
1202
C**   PARAMETER:   PHIS     REAL   GEOGR. BREITE DES PUNKTES IM ROT.SYS.
1203
C**                LAMS     REAL   GEOGR. LAENGE DES PUNKTES IM ROT.SYS.
1204
C**                POLPHI   REAL   WAHRE GEOGR. BREITE DES NORDPOLS
1205
C**                POLLAM   REAL   WAHRE GEOGR. LAENGE DES NORDPOLS
1206
C**   AUSGABE-
1207
C**   PARAMETER:   WAHRE GEOGRAPHISCHE BREITE ALS WERT DER FUNKTION
1208
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
1209
C**
1210
C**   COMMON-
1211
C**   BLOECKE  :   KEINE
1212
C**
1213
C**   FEHLERBE-
1214
C**   HANDLUNG :   KEINE
1215
C**   VERFASSER:   D.MAJEWSKI
1216
 
1217
      REAL        LAMS,PHIS,POLPHI,POLLAM
1218
 
1219
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
1220
 
1221
      SINPOL = SIN(ZPIR18*POLPHI)
1222
      COSPOL = COS(ZPIR18*POLPHI)
1223
      ZPHIS  = ZPIR18*PHIS
1224
      ZLAMS  = LAMS
1225
      IF(ZLAMS.GT.180.0) ZLAMS = ZLAMS - 360.0
1226
      ZLAMS  = ZPIR18*ZLAMS
1227
      ARG     = COSPOL*COS(ZPHIS)*COS(ZLAMS) + SINPOL*SIN(ZPHIS)
1228
 
1229
      PHSTOPH = ZRPI18*ASIN(ARG)
1230
 
1231
      RETURN
1232
      END
1233
 
1234
 
1235
      REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)
1236
C
1237
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
1238
C
1239
C**** LMTOLMS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM
1240
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
1241
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
1242
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1243
C**   AUFRUF   :   LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)
1244
C**   ENTRIES  :   KEINE
1245
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF
1246
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
1247
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1248
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1249
C**   VERSIONS-
1250
C**   DATUM    :   03.05.90
1251
C**
1252
C**   EXTERNALS:   KEINE
1253
C**   EINGABE-
1254
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
1255
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
1256
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
1257
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
1258
C**   AUSGABE-
1259
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
1260
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
1261
C**
1262
C**   COMMON-
1263
C**   BLOECKE  :   KEINE
1264
C**
1265
C**   FEHLERBE-
1266
C**   HANDLUNG :   KEINE
1267
C**   VERFASSER:   G. DE MORSIER
1268
 
1269
      REAL        LAM,PHI,POLPHI,POLLAM
1270
 
1271
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
1272
 
1273
      ZSINPOL = SIN(ZPIR18*POLPHI)
1274
      ZCOSPOL = COS(ZPIR18*POLPHI)
1275
      ZLAMPOL =     ZPIR18*POLLAM
1276
      ZPHI    =     ZPIR18*PHI
1277
      ZLAM    = LAM
1278
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
1279
      ZLAM    = ZPIR18*ZLAM
1280
 
1281
      ZARG1   = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI)
1282
      ZARG2   = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI)
1283
      IF (ABS(ZARG2).LT.1.E-30) THEN
1284
        IF (ABS(ZARG1).LT.1.E-30) THEN
1285
          LMTOLMS =   0.0
1286
        ELSEIF (ZARG1.GT.0.) THEN
1287
              LMTOLMS =  90.0
1288
            ELSE
1289
              LMTOLMS = -90.0
1290
            ENDIF
1291
      ELSE
1292
        LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
1293
      ENDIF
1294
 
1295
      RETURN
1296
      END
1297
 
1298
 
1299
      REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
1300
C
1301
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
1302
C
1303
C**** PHTOPHS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
1304
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
1305
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
1306
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1307
C**   AUFRUF   :   PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)
1308
C**   ENTRIES  :   KEINE
1309
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF
1310
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
1311
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
1312
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
1313
C**   VERSIONS-
1314
C**   DATUM    :   03.05.90
1315
C**
1316
C**   EXTERNALS:   KEINE
1317
C**   EINGABE-
1318
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
1319
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
1320
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
1321
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
1322
C**   AUSGABE-
1323
C**   PARAMETER:   ROTIERTE BREITE PHIS ALS WERT DER FUNKTION
1324
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
1325
C**
1326
C**   COMMON-
1327
C**   BLOECKE  :   KEINE
1328
C**
1329
C**   FEHLERBE-
1330
C**   HANDLUNG :   KEINE
1331
C**   VERFASSER:   G. DE MORSIER
1332
 
1333
      REAL        LAM,PHI,POLPHI,POLLAM
1334
 
1335
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
1336
 
1337
      ZSINPOL = SIN(ZPIR18*POLPHI)
1338
      ZCOSPOL = COS(ZPIR18*POLPHI)
1339
      ZLAMPOL = ZPIR18*POLLAM
1340
      ZPHI    = ZPIR18*PHI
1341
      ZLAM    = LAM
1342
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
1343
      ZLAM    = ZPIR18*ZLAM
1344
      ZARG    = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI)
1345
 
1346
      PHTOPHS = ZRPI18*ASIN(ZARG)
1347
 
1348
      RETURN
1349
      END
3 michaesp 1350