Subversion Repositories pvinversion.ecmwf

      PROGRAM coastline
      
c     *********************************************************************
c     * Add the coastline and land sea mask to the REF file               *
c     * Michael Sprenger / Winter 2006,07                                 *
c     *********************************************************************

      implicit none

c     --------------------------------------------------------------------
c     Declaration of parameters and variables
c     --------------------------------------------------------------------

c     Input parameters
      character*80 coastfile
      character*80 reffile
      real         clon,clat,crot
      
c     Output arrays
      integer      nmax
      parameter    (nmax=100000)
      real,allocatable, dimension (:,:)   :: landsea
      real,allocatable, dimension (:,:)   :: xcor,ycor
      real         xpos(nmax),ypos(nmax)
      real         lon(nmax), lat(nmax)
      real         rlon(nmax),rlat(nmax)
      integer      ncoast
      integer      rnx,rny
      real         rxmin,rymin,rxmax,rymax
      real         rdx,rdy

c     Auxiliary variables
      integer      i,j
      character*80 name
      integer      cdfid,cstid
      character*80 cfn
      integer      ierr,stat
      character*80 varname
      real         varmin(4),varmax(4),stag(4)
      integer      vardim(4),ndim  
      real         mdv
      integer      isok
      integer      nvars
      integer      ntimes
      real         time
      character*80 vnam(100)

c     --------------------------------------------------------------------
c     Preparations
c     --------------------------------------------------------------------

      print*,'*********************************************************'
      print*,'* coastline                                             *'
      print*,'*********************************************************'


c     Read parameter file
      open(10,file='fort.10')
       read(10,*) reffile
       read(10,*) coastfile
       read(10,*) name,clon
       if ( trim(name).ne.'CLON' ) stop
       read(10,*) name,clat
       if ( trim(name).ne.'CLAT' ) stop      
       read(10,*) name,crot
       if ( trim(name).ne.'CROT' ) stop    
      close(10)

      print*
      print*,trim(reffile)
      print*,trim(coastfile)
      print*,clon,clat,crot

c     Read grid parameters
      call cdfopn(reffile,cdfid,ierr)
      if (ierr.ne.0) goto 998
      call getvars(cdfid,nvars,vnam,ierr)
      if (ierr.ne.0) goto 998     
      varname='X'      
      call getdef(cdfid,varname,ndim,mdv,vardim,
     >            varmin,varmax,stag,ierr)
      if (ierr.ne.0) goto 998
      call gettimes(cdfid,time,ntimes,ierr)
      if (ierr.ne.0) goto 998
      call clscdf(cdfid,ierr)
      rnx=vardim(1)
      rny=vardim(2)
      rxmin=varmin(1)
      rymin=varmin(2)
      rxmax=varmax(1)
      rymax=varmax(2)
      rdx=(rxmax-rxmin)/real(rnx-1)
      rdy=(rymax-rymin)/real(rny-1)

c     Allocate memory
      allocate(landsea(rnx,rny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array landsea ***'
      allocate(xcor(rnx,rny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array xcor ***'
      allocate(ycor(rnx,rny),stat=stat)
      if (stat.ne.0) print*,'*** error allocating array ycor ***'

c     Read data
      call cdfopn(reffile,cdfid,ierr)
      if (ierr.ne.0) goto 998
      varname='X'
      call getdat(cdfid,varname,time,1,xcor,ierr)
      if (ierr.ne.0) goto 998
      varname='Y'
      call getdat(cdfid,varname,time,1,ycor,ierr)
      if (ierr.ne.0) goto 998
      call clscdf(cdfid,ierr)

c     Read the coastfile
      open(10,file=coastfile)
      
       ncoast=1
 100   read(10,*,end=200) lon(ncoast),lat(ncoast)
       ncoast=ncoast+1
       goto 100

 200  close(10)
      ncoast=ncoast-1

c     Handle missing data values
      mdv=-100000.
      do i=ncoast,nmax
         lon(i)=mdv
         lat(i)=mdv
      enddo

c     --------------------------------------------------------------------
c     Transform coast(lat/lon) to coast(x/y)
c     --------------------------------------------------------------------
      
      call  getenvir (clon,clat,crot,
     >                xcor,ycor,rnx,rny,rxmin,rymin,rdx,rdy,mdv,
     >                lon,lat,rlon,rlat,xpos,ypos,nmax)

c     --------------------------------------------------------------------
c     Write output 
c     --------------------------------------------------------------------

c     Open output file and set some general parameters
      call cdfwopn(reffile,cdfid,ierr)
      if (ierr.ne.0) goto 998
      vardim(1)=1
      vardim(2)=1
      vardim(3)=nmax
      ndim=3

c     Write LAT
      isok=0
      varname='COAST_LAT'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,ndim,mdv,vardim,
     >               varmin,varmax,stag,ierr)
         if (ierr.ne.0) goto 998
      endif
      call putdat(cdfid,varname,time,0,lat,ierr)
      if (ierr.ne.0) goto 998
      print*,'W ',trim(varname),' ',trim(reffile)

c     Write LON
      isok=0
      varname='COAST_LON'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,ndim,mdv,vardim,
     >               varmin,varmax,stag,ierr)
         if (ierr.ne.0) goto 998
      endif
      call putdat(cdfid,varname,time,0,lon,ierr)
      if (ierr.ne.0) goto 998
      print*,'W ',trim(varname),' ',trim(reffile)

c     Write RLAT
      isok=0
      varname='COAST_RLAT'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,ndim,mdv,vardim,
     >               varmin,varmax,stag,ierr)
         if (ierr.ne.0) goto 998
      endif
      call putdat(cdfid,varname,time,0,rlat,ierr)
      if (ierr.ne.0) goto 998
      print*,'W ',trim(varname),' ',trim(reffile)

c     Write RLON
      isok=0
      varname='COAST_RLON'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,ndim,mdv,vardim,
     >               varmin,varmax,stag,ierr)
         if (ierr.ne.0) goto 998
      endif
      call putdat(cdfid,varname,time,0,rlon,ierr)
      if (ierr.ne.0) goto 998
      print*,'W ',trim(varname),' ',trim(reffile)

c     Write Y
      isok=0
      varname='COAST_Y'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,ndim,mdv,vardim,
     >               varmin,varmax,stag,ierr)
         if (ierr.ne.0) goto 998
      endif
      call putdat(cdfid,varname,time,0,ypos,ierr)
      if (ierr.ne.0) goto 998
      print*,'W ',trim(varname),' ',trim(reffile)

c     Write X
      isok=0
      varname='COAST_X'
      call check_varok(isok,varname,vnam,nvars)
      if (isok.eq.0) then
         call putdef(cdfid,varname,ndim,mdv,vardim,
     >               varmin,varmax,stag,ierr)
         if (ierr.ne.0) goto 998
      endif
      call putdat(cdfid,varname,time,0,xpos,ierr)
      if (ierr.ne.0) goto 998
      print*,'W ',trim(varname),' ',trim(reffile)

c     Close output file
      call clscdf(cdfid,ierr)
      if (ierr.ne.0) goto 998


c     --------------------------------------------------------------------
c     Exception handling
c     --------------------------------------------------------------------
      
      stop

 998  print*,'Problems with input netcdf file'
      stop

      END

      
c     --------------------------------------------------------------------------------
c     Subroutine to get environment of strcof
c     --------------------------------------------------------------------------------

      SUBROUTINE getenvir (clon,clat,rotation,
     >                     xcor,ycor,rnx,rny,rxmin,rymin,rdx,rdy,mdv,
     >                     lon,lat,rlon,rlat,xpos,ypos,n)

c     Rotate from a local quasi-cartesian coordiante system into lat/lon coordinate 
c     system.

      implicit none

c     Declaration of input and output parameters
      integer     n
      real        clon,clat,rotation
      real        lon(n), lat(n)
      real        rlon(n),rlat(n)
      real        xpos(n),ypos(n)
      integer     rnx,rny
      real        xcor(rnx,rny),ycor(rnx,rny)
      real        rxmin,rymin,rdx,rdy
      real        mdv

c     Set numerical and physical constants
      real        g2r
      parameter   (g2r=0.0174533)
      real        pi180
      parameter   (pi180=3.14159265359/180.)
      real        eps
      parameter   (eps=0.0001)

c     Auxiliary variables 
      real         pollon,pollat
      integer      i,j
      real         xind,yind
      real         rindx,rindy
      integer      indx,indy,indr,indu
      real         frac0i,frac0j,frac1i,frac1j
      real         rlon1(n),rlat1(n)

c     Externals
      real     lmtolms,phtophs
      external lmtolms,phtophs

c     First rotation
      pollon=clon-180.
      if (pollon.lt.-180.) pollon=pollon+360.
      pollat=90.-clat
      do i=1,n

         if ( (abs(lon(i)-mdv).gt.eps).and.
     >        (abs(lat(i)-mdv).gt.eps) ) then

c           First rotation
            pollon=clon-180.
            if (pollon.lt.-180.) pollon=pollon+360.
            pollat=90.-clat
            rlon1(i)=lmtolms(lat(i),lon(i),pollat,pollon)
            rlat1(i)=phtophs(lat(i),lon(i),pollat,pollon)            

c           Second coordinate transformation
            pollon=-180.
            pollat=90.+rotation
            rlon(i)=90.+lmtolms(rlat1(i),rlon1(i)-90.,pollat,pollon)
            rlat(i)=phtophs(rlat1(i),rlon1(i)-90.,pollat,pollon)   

c           Get rotated latitude and longitude
 100        if (rlon(i).lt.rxmin) then
               rlon(i)=rlon(i)+360.
               goto 100
            endif
 102        if (rlon(i).gt.(rxmin+real(rnx-1)*rdx)) then
               rlon(i)=rlon(i)-360.
               goto 102
            endif
            if ( (rlon(i).lt.rxmin).or.
     >           (rlon(i).gt.(rxmin+real(rnx-1)*rdx)).or.
     >           (rlat(i).lt.rymin).or.
     >           (rlat(i).gt.(rymin+real(rny-1)*rdy)) ) then
               
               rlon(i)=mdv
               rlat(i)=mdv
            endif

c           Get x and y coordinates
            rindx=(rlon(i)-rxmin)/rdx+1.
            rindy=(rlat(i)-rymin)/rdy+1.
            indx=int(rindx)
            indy=int(rindy)
            if ((indx.gt.1).and.(indx.lt.rnx).and.
     >          (indy.gt.1).and.(indy.lt.rny)) then

               indr=indx+1
               if (indr.gt.rnx) indr=1
               indu=indy+1
               if (indu.gt.rny) indu=rny

               frac0i=rindx-float(indx)
               frac0j=rindy-float(indy)
               frac1i=1.-frac0i
               frac1j=1.-frac0j
               xpos(i) = xcor(indx ,indy) * frac1i * frac1j
     &                 + xcor(indx ,indu) * frac1i * frac0j
     &                 + xcor(indr ,indy) * frac0i * frac1j
     &                 + xcor(indr ,indu) * frac0i * frac0j
               ypos(i) = ycor(indx ,indy) * frac1i * frac1j
     &                 + ycor(indx ,indu) * frac1i * frac0j
     &                 + ycor(indr ,indy) * frac0i * frac1j
     &                 + ycor(indr ,indu) * frac0i * frac0j

            else

               xpos(i)=mdv
               ypos(i)=mdv

            endif

         else
            
            rlon(i)=mdv
            rlat(i)=mdv
            xpos(i)=mdv
            ypos(i)=mdv
            
         endif

      enddo

      END    

      
c     --------------------------------------------------------------------------------
c     Transformation routine: LMSTOLM and PHSTOPH from library gm2em
c     --------------------------------------------------------------------------------

      REAL FUNCTION LMTOLMS (PHI, LAM, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
C
C**** LMTOLMS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   LAM = LMTOLMS (PHI, LAM, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE LAM AUF
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
C**   AUSGABE-
C**   PARAMETER:   WAHRE GEOGRAPHISCHE LAENGE ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   G. DE MORSIER
 
      REAL        LAM,PHI,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL =     ZPIR18*POLLAM
      ZPHI    =     ZPIR18*PHI
      ZLAM    = LAM
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
      ZLAM    = ZPIR18*ZLAM
 
      ZARG1   = - SIN(ZLAM-ZLAMPOL)*COS(ZPHI)
      ZARG2   = - ZSINPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL)+ZCOSPOL*SIN(ZPHI)
      IF (ABS(ZARG2).LT.1.E-30) THEN
        IF (ABS(ZARG1).LT.1.E-30) THEN
          LMTOLMS =   0.0
        ELSEIF (ZARG1.GT.0.) THEN
              LMTOLMS =  90.0
            ELSE
              LMTOLMS = -90.0
            ENDIF
      ELSE
        LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
      ENDIF
 
      RETURN
      END


      REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
C
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
C
C**** PHTOPHS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
C****                 AUF EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS)
C****                 IM ROTIERTEN SYSTEM. DER NORDPOL DES SYSTEMS HAT
C****                 DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   AUFRUF   :   PHI = PHTOPHS (PHI, LAM, POLPHI, POLLAM)
C**   ENTRIES  :   KEINE
C**   ZWECK    :   UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI AUF
C**                EINEM PUNKT MIT DEN KOORDINATEN (PHIS, LAMS) IM
C**                ROTIERTEN SYSTEM. DER NORDPOL DIESES SYSTEMS HAT
C**                DIE WAHREN KOORDINATEN (POLPHI, POLLAM)
C**   VERSIONS-
C**   DATUM    :   03.05.90
C**
C**   EXTERNALS:   KEINE
C**   EINGABE-
C**   PARAMETER:   PHI    REAL BREITE DES PUNKTES IM GEOGR. SYSTEM
C**                LAM    REAL LAENGE DES PUNKTES IM GEOGR. SYSTEM
C**                POLPHI REAL GEOGR.BREITE DES N-POLS DES ROT. SYSTEMS
C**                POLLAM REAL GEOGR.LAENGE DES N-POLS DES ROT. SYSTEMS
C**   AUSGABE-
C**   PARAMETER:   ROTIERTE BREITE PHIS ALS WERT DER FUNKTION
C**                ALLE WINKEL IN GRAD (NORDEN>0, OSTEN>0)
C**
C**   COMMON-
C**   BLOECKE  :   KEINE
C**
C**   FEHLERBE-
C**   HANDLUNG :   KEINE
C**   VERFASSER:   G. DE MORSIER
 
      REAL        LAM,PHI,POLPHI,POLLAM
 
      DATA        ZRPI18 , ZPIR18  / 57.2957795 , 0.0174532925 /
 
      ZSINPOL = SIN(ZPIR18*POLPHI)
      ZCOSPOL = COS(ZPIR18*POLPHI)
      ZLAMPOL = ZPIR18*POLLAM
      ZPHI    = ZPIR18*PHI
      ZLAM    = LAM
      IF(ZLAM.GT.180.0) ZLAM = ZLAM - 360.0
      ZLAM    = ZPIR18*ZLAM
      ZARG    = ZCOSPOL*COS(ZPHI)*COS(ZLAM-ZLAMPOL) + ZSINPOL*SIN(ZPHI)
 
      PHTOPHS = ZRPI18*ASIN(ZARG)
 
      RETURN
      END

     
c     ----------------------------------------------------------------
c     Check whether variable is found on netcdf file
c     ----------------------------------------------------------------

      subroutine check_varok (isok,varname,varlist,nvars)

c     Check whether the variable <varname> is in the list <varlist(nvars)>.
c     If this is the case, <isok> is incremented by 1. Otherwise <isok>
c     keeps its value.

      implicit none

c     Declaraion of subroutine parameters
      integer      isok
      integer      nvars
      character*80 varname
      character*80 varlist(nvars)

c     Auxiliary variables
      integer      i

c     Main
      do i=1,nvars
         if (trim(varname).eq.trim(varlist(i))) isok=isok+1
      enddo

      end