Subversion Repositories lagranto.ecmwf

Rev

Rev 3 | Rev 9 | Go to most recent revision | Show entire file | Regard whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 3 Rev 5
Line 68... Line 68...
68 
      character*80                           timecheck       ! Either 'yes' or 'no'
 68 
      character*80                           timecheck       ! Either 'yes' or 'no'
69 
      character*80                           intmode         ! Interpolation mode ('normal', 'nearest')
 69 
      character*80                           intmode         ! Interpolation mode ('normal', 'nearest')
70 
	  real                                   pmin,pmax       ! Pressure range for output grid
 70 
	  real                                   pmin,pmax       ! Pressure range for output grid
71 
	  integer                                npre            ! Number of pressure levels in output grid
 71 
	  integer                                npre            ! Number of pressure levels in output grid
72 
	  character*80                           centering       ! Centering around trajectory position ('yes','no')
 72 
	  character*80                           centering       ! Centering around trajectory position ('yes','no')
- 
 
 73 
	  character*80                       direction       ! Direction of lidar (vertical,lat,lon,normal)
- 
 
 74 
      character*80                           dumpcoord       ! Dumping coordinates ('yes','no')
73 
 
 75 
 
74 
c     Trajectories
 76 
c     Trajectories
75 
      real,allocatable, dimension (:,:,:) :: trainp          ! Input trajectories (ntra,ntim,ncol)
 77 
      real,allocatable, dimension (:,:,:) :: trainp                ! Input trajectories (ntra,ntim,ncol)
76 
      integer                                reftime(6)      ! Reference date
 78 
      integer                                reftime(6)            ! Reference date
77 
      character*80                           varsinp(100)    ! Field names for input trajectory
 79 
      character*80                           varsinp(100)          ! Field names for input trajectory
78 
      integer                                fid,fod         ! File identifier for inp and out trajectories
 80 
      integer                                fid,fod               ! File identifier for inp and out trajectories
79 
      real                                   x0,y0,p0        ! Position of air parcel (physical space)
 81 
      real                                   x0_tra,y0_tra,p0_tra  ! Position of air parcel (physical space)
80 
      real                                   reltpos0        ! Relative time of air parcel
 82 
      real                                   reltpos0              ! Relative time of air parcel
81 
      real                                   xind,yind,pind  ! Position of air parcel (grid space)
 83 
      real                                   xind,yind,pind        ! Position of air parcel (grid space)
82 
      integer                                fbflag          ! Flag for forward (1) or backward (-1) trajectories
 84 
      integer                                fbflag                ! Flag for forward (1) or backward (-1) trajectories
83 
 
 85 
 
84 
c     Meteorological fields from input file
 86 
c     Meteorological fields from input file
Line 137... Line 139...
137 
      real                                   time
 139 
      real                                   time
138 
      character*80                           longname
 140 
      character*80                           longname
139 
      character*80                           unit
 141 
      character*80                           unit
140 
      integer                                ind_time
 142 
      integer                                ind_time
141 
      integer                                ind_pre
 143 
      integer                                ind_pre
- 
 
 144 
      real                                   rlat,rlon
- 
 
 145 
      real                                   x0,y0,p0
- 
 
 146 
      real                                   vx0,vy0,vx1,vy1
- 
 
 147 
      real                                   rotation,lon,lat
142 
 
 148 
 
143 
c     Externals
149 
c     Externals
144 
      real                                   int_index4
 150 
      real                                   int_index4
145 
      external                               int_index4
 151 
      external                               int_index4
146 
 
 152 
 
Line 191... Line 197...
191 
       enddo
 197 
       enddo
192 
       read(10,*) timecheck
 198 
       read(10,*) timecheck
193 
       read(10,*) intmode
 199 
       read(10,*) intmode
194 
       read(10,*) pmin,pmax,npre
 200 
       read(10,*) pmin,pmax,npre
195 
       read(10,*) centering
 201 
       read(10,*) centering
- 
 
 202 
       read(10,*) direction
- 
 
 203 
       read(10,*) dumpcoord
196 
      close(10)
 204 
      close(10)
197 
 
 205 
 
- 
 
 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 
 
198 
c     Remove commented tracing fields
 216 
c     Remove commented tracing fields
199 
      itrace0 = 1
 217 
      itrace0 = 1
200 
      do while ( itrace0.le.ntrace0)
218 
      do while ( itrace0.le.ntrace0)
201 
         string = tvar(itrace0)
 219 
         string = tvar(itrace0)
202 
         if ( string(1:1).eq.'#' ) then
 220 
         if ( string(1:1).eq.'#' ) then
Line 309... Line 327...
309 
     >                trim(tvar(i)),' : online calc not supported'
 327 
     >                trim(tvar(i)),' : online calc not supported'
310 
         endif
 328 
         endif
311 
      enddo
 329 
      enddo
312 
      print*,'  Output (pmin,pmax,n)   : ',pmin,pmax,npre
 330 
      print*,'  Output (pmin,pmax,n)   : ',pmin,pmax,npre
313 
      print*,'  Centering              : ',trim(centering)
 331 
      print*,'  Centering              : ',trim(centering)
- 
 
 332 
      print*,'  Orientation            : ',trim(direction)
- 
 
 333 
      print*,'  Coordinate Dump        : ',trim(dumpcoord)
314 
      print*
 334 
      print*
315 
      print*,'---- INPUT DATA FILES -----------------------------------'
 335 
      print*,'---- INPUT DATA FILES -----------------------------------'
316 
      print*
 336 
      print*
317 
      call frac2hhmm(tstart,tload)
 337 
      call frac2hhmm(tstart,tload)
318 
      print*,'  Time of 1st data file  : ',tload
 338 
      print*,'  Time of 1st data file  : ',tload
Line 394... Line 414...
394 
         print*,'     correspond to the reference date. Otherwise'
 414 
         print*,'     correspond to the reference date. Otherwise'
395 
         print*,'     the tracing will give wrong results... STOP'
 415 
         print*,'     the tracing will give wrong results... STOP'
396 
         stop
 416 
         stop
397 
      endif
 417 
      endif
398
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 
 
399 
c     --------------------------------------------------------------------
 424 
c     --------------------------------------------------------------------
400 
c     Trace the fields (fields available on input files)
 425 
c     Trace the fields (fields available on input files)
401 
c     --------------------------------------------------------------------
 426 
c     --------------------------------------------------------------------
402 
 
 427 
 
403 
      print*
 428 
      print*
Line 510... Line 535...
510 
             endif
 535 
             endif
511 
 
 536 
 
512 
c            Loop over all trajectories
 537 
c            Loop over all trajectories
513 
             do k=1,ntra
 538 
             do k=1,ntra
514
539
515 
c               Set the horizontal position where to interpolate to
 540 
c               Set the trajectory position
516 
                x0       = trainp(k,j,2)                          ! Longitude
 541 
                x0_tra = trainp(k,j,2)                          ! Longitude
517 
                y0       = trainp(k,j,3)                          ! Latitude
 542 
                y0_tra = trainp(k,j,3)                          ! Latitude
- 
 
 543 
                p0_tra = trainp(k,j,4)                          ! Pressure
- 
 
 544 
 
- 
 
 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
518 
 
 572 
 
519 
c               Set the relative time
 573 
c               Set the relative time
520 
                call hhmm2frac(trainp(k,j,1),tfrac)
 574 
                call hhmm2frac(trainp(k,j,1),tfrac)
521 
                reltpos0 = fbflag * (tfrac-time0)/timeinc
575 
                reltpos0 = fbflag * (tfrac-time0)/timeinc
522
576 
 
- 
 
 577 
c               Loop over pressure profile (or other positions for horizontal mode)
- 
 
 578 
	        do l=1,npre
- 
 
 579 
 
- 
 
 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
- 
 
 588 
 
- 
 
 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
- 
 
 606 
 
- 
 
 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 
 
523 
c               Handle periodic boundaries in zonal direction
 624 
c                     Handle periodic boundaries in zonal direction
524 
                if ( (x0.gt.xmax).and.(per.ne.0) ) x0 = x0 - 360.
 625 
                      if ( (x0.gt.xmax).and.(per.ne.0) ) x0 = x0 - 360.
525 
                if ( (x0.lt.xmin).and.(per.ne.0) ) x0 = x0 + 360.
 626 
                      if ( (x0.lt.xmin).and.(per.ne.0) ) x0 = x0 + 360.
526 
 
 627 
 
527 
c               Handle pole problems for hemispheric data (taken from caltra.f)
 628 
c                     Handle pole problems for hemispheric data (taken from caltra.f)
Line 535... Line 636...
535 
                endif
 636 
                      endif
536 
                if (y0.gt.89.99) then
 637 
                      if (y0.gt.89.99) then
537 
                   y0=89.99
 638 
                         y0=89.99
538 
                endif
 639 
                      endif
539 
 
 640 
 
540 
c               Loop over pressure profile
 - 
 
541 
	            do l=1,npre
 - 
 
542 
 
 - 
 
543 
c                 Set the pressure where to interpolate to
 641 
c                 If requested, dump the lidar coordinates
544 
	              p0 = pmin + real(l-1)/real(npre-1) * (pmax-pmin)
 642 
                  if ( (dumpcoord.eq.'yes').and.(i.eq.1) ) then
545 
	              if ( centering.eq.'yes' )then
 643 
                     write(10,'3f10.2') x0,y0,trainp(k,j,1)
546 
	              	  p0 = p0 + trainp(k,j,4)
 644 
                     write(10,'3f10.2') x0_tra,y0_tra,5.
547 
	              endif
 645 
                  endif
548 
 
 646 
 
549 
C                 Get the index where to interpolate (x0,y0,p0)
 647 
C                 Get the index where to interpolate (x0,y0,p0)
550 
                  if ( (abs(x0-mdv).gt.eps).and.
 648 
                  if ( (abs(x0-mdv).gt.eps).and.
551 
     >                 (abs(y0-mdv).gt.eps) )
 649 
     >                 (abs(y0-mdv).gt.eps) )
Line 597... Line 695...
597 
	              endif
 695 
	              endif
598 
 
 696 
 
599 
c              End loop over all pressure levels
 697 
c              End loop over all pressure levels
600 
	           enddo
 698 
	           enddo
601 
 
 699 
 
602 
c	           Save the output trajectory position
 700 
c	           Save output - time index
603 
	           ind_time = j
 701 
	           ind_time = j
- 
 
 702 
 
- 
 
 703 
c                  Save output - space index for 'no centering'
604 
	           if ( centering.eq.'no' ) then
 704 
	           if ( centering.eq.'no' ) then
- 
 
 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
605 
	              ind_pre  = nint( real(npre) *
 709 
                         ind_pre  = nint( real(npre) *
606 
     >          	       ( (trainp(k,j,4) - pmin)/(pmax-pmin) ) + 1.)
 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'
607 
	           else
 717 
	           else
608 
	           	  ind_pre  = nint( real(npre) *
 718 
	           	  ind_pre  = nint( real(npre) *
609 
     >          	       ( (0.            - pmin)/(pmax-pmin) ) + 1.)
 719 
     >          	       ( (0.            - pmin)/(pmax-pmin) ) + 1.)
610 
	           endif
 720 
	           endif
611 
 
 721 
 
- 
 
 722 
c                  Update the output array
612 
	           if ( (ind_time.ge.1).and.(ind_time.le.ntim).and.
 723 
	           if ( (ind_time.ge.1).and.(ind_time.le.ntim).and.
613 
     >              (ind_pre .ge.1).and.(ind_pre .le.npre) )
 724 
     >              (ind_pre .ge.1).and.(ind_pre .le.npre) )
614 
     >         then
 725 
     >         then
615 
                    out_pos(ind_time,ind_pre) =
 726 
                    out_pos(ind_time,ind_pre) =
616 
     >                          	out_pos(ind_time,ind_pre) + 1.
 727 
     >                          	out_pos(ind_time,ind_pre) + 1.
Line 635... Line 746...
635 
              do k=1,ntim
 746 
              do k=1,ntim
636 
                 times(k) = trainp(1,k,1)
 747 
                 times(k) = trainp(1,k,1)
637 
              enddo
 748 
              enddo
638 
              call writecdf2D_cf
 749 
              call writecdf2D_cf
639 
     >            (cdfname,varname,longname,unit,out_pos,time,levels,
 750 
     >            (cdfname,varname,longname,unit,out_pos,time,levels,
640 
     >             times,npre,ntim,1,1)
 751 
     >             times,npre,ntim,1,1,direction)
641 
	      endif
 752 
	      endif
642 
 
 753 
 
643 
c	      If no valid lidar count: set the field to missing data
 754 
c	      If no valid lidar count: set the field to missing data
644 
          do k=1,ntim
 755 
          do k=1,ntim
645 
          	do l=1,npre
 756 
          	do l=1,npre
Line 672... Line 783...
672 
	      longname = 'sum over all '//trim(tvar(i))//' profiles'
 783 
	      longname = 'sum over all '//trim(tvar(i))//' profiles'
673 
	      unit     = 'not given'
 784 
	      unit     = 'not given'
674 
	      time     = 0.
 785 
	      time     = 0.
675 
          call writecdf2D_cf
 786 
          call writecdf2D_cf
676 
     >            (cdfname,varname,longname,unit,out_val,time,levels,
 787 
     >            (cdfname,varname,longname,unit,out_val,time,levels,
677 
     >             times,npre,ntim,0,1)
 788 
     >             times,npre,ntim,0,1,direction)
678 
 
 789 
 
679 
	  	  cdfname  = outfile
 790 
	  	  cdfname  = outfile
680 
	      varname  = trim(tvar(i))//'_CNT'
 791 
	      varname  = trim(tvar(i))//'_CNT'
681 
	      longname = 'counts of all '//trim(tvar(i))//' profiles'
 792 
	      longname = 'counts of all '//trim(tvar(i))//' profiles'
682 
	      unit     = 'not given'
 793 
	      unit     = 'not given'
683 
	      time     = 0.
 794 
	      time     = 0.
684 
          call writecdf2D_cf
 795 
          call writecdf2D_cf
685 
     >            (cdfname,varname,longname,unit,out_cnt,time,levels,
 796 
     >            (cdfname,varname,longname,unit,out_cnt,time,levels,
686 
     >             times,npre,ntim,0,1)
 797 
     >             times,npre,ntim,0,1,direction)
687 
 
 798 
 
688 
c         Exit point for loop over all tracing variables
 799 
c         Exit point for loop over all tracing variables
689 
 110      continue
 800 
 110      continue
690 
 
 801 
 
691 
c	   End loop over all lidar variables
 802 
c	   End loop over all lidar variables
Line 699... Line 810...
699 
c     Write status information
 810 
c     Write status information
700 
      print*
 811 
      print*
701 
      print*,'---- WRITE OUTPUT LIDAR FIELDS --------------------------'
 812 
      print*,'---- WRITE OUTPUT LIDAR FIELDS --------------------------'
702 
      print*
 813 
      print*
703 
 
 814 
 
- 
 
 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
704 
 
 820 
 
705 
c     Write some status information, and end of program message
 821 
c     Write some status information, and end of program message
706 
      print*
822 
      print*
707 
      print*,'---- STATUS INFORMATION --------------------------------'
 823 
      print*,'---- STATUS INFORMATION --------------------------------'
708 
      print*
 824 
      print*
Line 723... Line 839...
723 
c     Subroutines to write 2D CF netcdf output file
 839 
c     Subroutines to write 2D CF netcdf output file
724 
c     --------------------------------------------------------------------
 840 
c     --------------------------------------------------------------------
725 
 
 841 
 
726 
      subroutine writecdf2D_cf
 842 
      subroutine writecdf2D_cf
727 
     >          (cdfname,varname,longname,unit,arr,time,levels,times,
 843 
     >          (cdfname,varname,longname,unit,arr,time,levels,times,
728 
     >           npre,ntim,crefile,crevar)
 844 
     >           npre,ntim,crefile,crevar,direction)
729 
 
 845 
 
730 
c     Create and write to the CF netcdf file <cdfname>. The variable
 846 
c     Create and write to the CF netcdf file <cdfname>. The variable
731 
c     with name <varname> and with time <time> is written. The data
 847 
c     with name <varname> and with time <time> is written. The data
732 
c     are in the two-dimensional array <arr>.  The flags <crefile> and
 848 
c     are in the two-dimensional array <arr>.  The flags <crefile> and
733 
c     <crevar> determine whether the file and/or the variable should
 849 
c     <crevar> determine whether the file and/or the variable should
Line 744... Line 860...
744 
      real         arr(ntim,npre)
 860 
      real         arr(ntim,npre)
745 
      real         levels(npre)
 861 
      real         levels(npre)
746 
      real         times (ntim)
 862 
      real         times (ntim)
747 
      real         time
 863 
      real         time
748 
      integer      crefile,crevar
 864 
      integer      crefile,crevar
- 
 
 865 
      character*80 direction
749 
 
 866 
 
750 
c     Numerical epsilon
 867 
c     Numerical epsilon
751 
      real         eps
 868 
      real         eps
752 
      parameter    (eps=1.e-5)
 869 
      parameter    (eps=1.e-5)
753 
 
 870 
 
Line 779... Line 896...
779 
 
 896 
 
780 
c     Define dimensions
 897 
c     Define dimensions
781 
      ierr=nf90_def_dim(ncID,'level',npre, LevDimID )
 898 
      ierr=nf90_def_dim(ncID,'level',npre, LevDimID )
782 
      ierr=nf90_def_dim(ncID,'time' ,ntim, TimeDimID)
 899 
      ierr=nf90_def_dim(ncID,'time' ,ntim, TimeDimID)
783 
 
 900 
 
784 
c     Define coordinate Variables
 901 
c     Define space coordinate
785 
      ierr = nf90_def_var(ncID,'level',NF90_FLOAT,
 902 
      ierr = nf90_def_var(ncID,'level',NF90_FLOAT,
786 
     >     (/ LevDimID /),varLevID)
 903 
     >     (/ LevDimID /),varLevID)
- 
 
 904 
      if ( direction.eq.'vertical' ) then
787 
      ierr = nf90_put_att(ncID, varLevID, "standard_name","level")
 905 
        ierr = nf90_put_att(ncID, varLevID, "standard_name","level")
788 
      ierr = nf90_put_att(ncID, varLevID, "units"        ,"hPa")
 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
789 
 
 917 
 
- 
 
 918 
c     Define time coordinate
790 
      ierr = nf90_def_var(ncID,'time',NF90_FLOAT,
 919 
      ierr = nf90_def_var(ncID,'time',NF90_FLOAT,
791 
     >     (/ TimeDimID /), varTimeID)
 920 
     >     (/ TimeDimID /), varTimeID)
792 
      ierr = nf90_put_att(ncID, varTimeID, "long_name",    "time")
 921 
      ierr = nf90_put_att(ncID, varTimeID, "long_name",    "time")
793 
      ierr = nf90_put_att(ncID, varTimeID, "units",       "hours")
 922 
      ierr = nf90_put_att(ncID, varTimeID, "units",       "hours")
794 
 
 923 
 
Line 820... Line 949...
820 
 
 949 
 
821 
c     ---- Define a new variable - skip if <crevar=0> -----------------------
 950 
c     ---- Define a new variable - skip if <crevar=0> -----------------------
822 
 
 951 
 
823 
      if ( crevar.ne.1 ) goto 110
 952 
      if ( crevar.ne.1 ) goto 110
824 
 
 953 
 
- 
 
 954 
      print*,'Now defining ',trim(varname)
- 
 
 955 
 
825 
c     Open the file for read/write access
 956 
c     Open the file for read/write access
826 
      ierr = nf90_open  (trim(cdfname), NF90_WRITE  , ncID)
 957 
      ierr = nf90_open  (trim(cdfname), NF90_WRITE  , ncID)
827 
 
 958 
 
828 
c     Get the IDs for dimensions
 959 
c     Get the IDs for dimensions
829 
      ierr = nf90_inq_dimid(ncID,'level', LevDimID )
 960 
      ierr = nf90_inq_dimid(ncID,'level', LevDimID )
Line 832... Line 963...
832 
c     Enter define mode
 963 
c     Enter define mode
833 
      ierr = nf90_redef(ncID)
 964 
      ierr = nf90_redef(ncID)
834 
 
 965 
 
835 
c     Write definition and add attributes
 966 
c     Write definition and add attributes
836 
      ierr = nf90_def_var(ncID,varname,NF90_FLOAT,
 967 
      ierr = nf90_def_var(ncID,varname,NF90_FLOAT,
837 
     >                   (/ varTimeID, varLevID /),varID)
 968 
     >                   (/ TimeDimID, LevDimID /),varID)
838 
      ierr = nf90_put_att(ncID, varID, "long_name" , longname )
 969 
      ierr = nf90_put_att(ncID, varID, "long_name" , longname )
839 
      ierr = nf90_put_att(ncID, varID, "units"     , unit     )
 970 
      ierr = nf90_put_att(ncID, varID, "units"     , unit     )
840 
      ierr = nf90_put_att(ncID, varID, '_FillValue', -999.99  )
 971 
      ierr = nf90_put_att(ncID, varID, '_FillValue', -999.99  )
841 
 
 972 
 
842 
c     Check whether definition was successful
 973 
c     Check whether definition was successful
Line 879... Line 1010...
879 
         write(*,*) 'in clscdf_CF'
 1010 
         write(*,*) 'in clscdf_CF'
880 
      endif
 1011 
      endif
881 
 
 1012 
 
882 
      end
 1013 
      end
883 
 
 1014 
 
- 
 
 1015 
c     ********************************************************************************
- 
 
 1016 
c     * Coordinate rotation - lidar normal to trajectory                             *
- 
 
 1017 
c     ********************************************************************************
- 
 
 1018 
 
- 
 
 1019 
c     --------------------------------------------------------------------------------
- 
 
 1020 
c     Backward coordinate transformation (Rotated lon/lat -> True lon/lat)
- 
 
 1021 
c     --------------------------------------------------------------------------------
- 
 
 1022 
 
- 
 
 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)
884 
 
 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
885 
 
 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
- 
 
 1092
- 
 
 1093 
         val1=vx1*vx2+vy1*vy2
- 
 
 1094 
         val2=-vy1*vx2+vx1*vy2
- 
 
 1095
- 
 
 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
886
1110
- 
 
 1111 
      angle=180./pi*val3
887
1112 
 
- 
 
 1113 
      END
888
1114 
 
- 
 
 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)
889 
 
 1345
- 
 
 1346 
      PHTOPHS = ZRPI18*ASIN(ZARG)
890 
 
 1347
- 
 
 1348 
      RETURN
- 
 
 1349 
      END
891
1350