WebSVN – lagranto.ecmwf – Diff – /trunk/create_startf/create_startf.f

 c     * whether to start from an isentropic or an isobaric         *
 c     * surface. The starting points are equidistantly distributed *
 c     * Michael Sprenger / Autumn 2004                             *
 c     **************************************************************
-      implicit none
+      use netcdf
+      implicit none
 c     --------------------------------------------------------------
 c     Set parameters
 c     --------------------------------------------------------------
       real              tfrac,frac
       real              radius,dist
       character*80      string
       character*80      selectstr
       character*80      umode_save
+      character*80      maskname
+      real              maskvalue
+      character*2       maskoper
+      integer           cdfid,varid,ierr
+      integer           nmask
+      integer           indx,indy
+      character         ch1,ch2
+      integer           len
       real,allocatable, dimension (:,:,:) :: fld0
       real,allocatable, dimension (:,:,:) :: fld1
       real,allocatable, dimension (:,:  ) :: sfc0
       real,allocatable, dimension (:,:)   :: sfc1
+      real,allocatable, dimension (:,:)   :: mask
 c     Externals
       real              int_index3     ! 3d interpolation
       external          int_index3
       real              sdis           ! Speherical distance
           read(10,*) lon1,lat1,radius
        elseif ( hmode.eq.'region.eqd' ) then        ! region: 2d equidistant
           read(10,*) iregion,ds
        elseif ( hmode.eq.'region.grid' ) then       ! iregion: 2d grid
           read(10,*) iregion
+       elseif ( hmode.eq.'mask.grid' ) then         ! 0/1 masked - grid points
+          read(10,*) hfile,maskname
+       elseif ( hmode.eq.'mask.eqd' ) then         ! 0/1 masked - euqidistant points
+          read(10,*) hfile,maskname,ds
        else
           print*,' ERROR: horizontal mode not supported ',trim(hmode)
           stop
        endif
          tfrac=tfrac/timeinc
       else
          tfrac=0.
       endif
+c     If a mask file is provided, no time interpolation is allowed
+      if ( (hmode.eq.'mask.grid').or.(hmode.eq.'mask.eqd') ) then
+         if ( abs(tfrac).gt.eps ) then
+            print*,' ERROR: no intermediate times allowed for ',
+     >             trim(hmode)
+            stop
+         endif
+      endif
 c     Read the region coordinates if needed
       if ( (hmode.eq.'region.eqd' ).or.
      >     (hmode.eq.'region.grid') ) then
          open(10,file=regionf)
      >         '      iregion                 : ',iregion
         write(*,'(a30,4f10.2)')
      >         '      xcorner                 : ',(xcorner(i),i=1,4)
         write(*,'(a30,4f10.2)')
      >         '      ycorner                 : ',(ycorner(i),i=1,4)
+      elseif ( hmode.eq.'mask.eqd' ) then
+        write(*,'(a30,a15,a15,f10.2)')
+     >         '      hfile,variable,ds       : ',trim(hfile),
+     >                                            trim(maskname),
+     >                                            ds
+      elseif ( hmode.eq.'mask.grid' ) then
+        write(*,'(a30,a15,a15)')
+     >         '      hfile,variable          : ',trim(hfile),
+     >                                            trim(maskname)
       endif
       print*
       print*,'  vmode                      : ',trim(vmode)
       if ( vmode.eq.'file') then
          print*,'      filename               : ',trim(vfile)
       elseif ( vmode.eq.'level' ) then
       elseif ( vmode.eq.'profile') then
          print*,'      lev1,lev2,n            : ',lev1,lev2,vn
       elseif ( vmode.eq.'grid') then
          print*,'      lev1,lev2              : ',lev1,lev2
       endif
       print*
       print*,'  umode                      : ',trim(umode)
       print*
       print*,'  time check                 : ',trim(timecheck)
       print*
+c     <mask.grid> and <mask.eqd>: split variable into name, operator,value
+      if ( (hmode.eq.'mask.grid').or.(hmode.eq.'mask.eqd') ) then
+          len = len_trim(maskname)
+          do i=1,len_trim(maskname)-1
+            ch1 = maskname(i:i)
+            ch2 = maskname(i+1:i+1)
+            if ( (ch1.eq.'<').and.(ch2.eq.'>') ) then
+               read(maskname(i+2:len),*) maskvalue
+               maskname = maskname(1:i-1)
+               maskoper = 'ne'
+               goto 90
+            elseif ( (ch1.eq.'<').and.(ch2.eq.'=') ) then
+               read(maskname(i+2:len),*) maskvalue
+               maskname = maskname(1:i-1)
+               maskoper = 'le'
+               goto 90
+            elseif ( (ch1.eq.'>').and.(ch2.eq.'=') ) then
+               read(maskname(i+2:len),*) maskvalue
+               maskname = maskname(1:i-1)
+               maskoper = 'ge'
+               goto 90
+            elseif ( (ch1.eq.'=').and.(ch2.eq.'=') ) then
+               read(maskname(i+2:len),*) maskvalue
+               maskname = maskname(1:i-1)
+               maskoper = 'eq'
+               goto 90
+            elseif ( ch1.eq.'=' ) then
+               read(maskname(i+1:len),*) maskvalue
+               maskname = maskname(1:i-1)
+               maskoper = 'eq'
+               goto 90
+            elseif ( ch1.eq.'>' ) then
+               read(maskname(i+1:len),*) maskvalue
+               maskname = maskname(1:i-1)
+               maskoper = 'gt'
+               goto 90
+            elseif ( ch1.eq.'<' ) then
+               read(maskname(i+1:len),*) maskvalue
+               maskname = maskname(1:i-1)
+               maskoper = 'lt'
+               goto 90
+             endif
+          enddo
+      endif
+continue
 c     ------------------------------------------------------------------
 c     Read grid parameters from inital files
 c     ------------------------------------------------------------------
 c     Get the time of the first and second data file
       if (stat.ne.0) print*,'*** error allocating array pv ***'
       allocate(pvs(nx,ny),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array pvs **'
       allocate(in(nx,ny,nz),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array in ***'
+      allocate(mask(nx,ny),stat=stat)
+      if (stat.ne.0) print*,'*** error allocating array mask ***'
 c     Allocate memory for temporary arrays for time interpolation
       allocate(fld0(nx,ny,nz),stat=stat)
       if (stat.ne.0) print*,'*** error allocating array tmp0 ***'
       allocate(fld1(nx,ny,nz),stat=stat)
                 pvs(i,j)=pv(i,j,1)
              enddo
           enddo
        endif
+c     ------ Load 0/1 label (mask) ------------------------------------
+      if ( (hmode.eq.'mask.grid').or.(hmode.eq.'mask.eqd') ) then
+c         Explicit netCDF calls - no consistency check for grids!
+          ierr = NF90_OPEN(hfile,nf90_nowrite, cdfid)
+          IF ( ierr /= nf90_NoErr ) PRINT *,NF90_STRERROR(ierr)
+          ierr = NF90_INQ_VARID(cdfid,maskname,varid)
+          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
+          ierr = nf90_get_var(cdfid,varid,mask)
+          IF(ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
+          ierr = NF90_CLOSE(cdfid)
+          IF( ierr /= nf90_NoErr) PRINT *,NF90_STRERROR(ierr)
+          nmask = 0
+          do i=1,nx
+            do j=1,ny
+                if ( maskoper.eq.'eq' ) then
+                   if ( abs(mask(i,j)-maskvalue).lt.eps ) then
+                      nmask = nmask + 1
+                      mask(i,j) = 1.
+                   else
+                      mask(i,j) = 0.
+                   endif
+                endif
+                if ( maskoper.eq.'ne' ) then
+                   if ( abs(mask(i,j)-maskvalue).gt.eps ) then
+                      nmask = nmask + 1
+                      mask(i,j) = 1.
+                   else
+                      mask(i,j) = 0.
+                   endif
+                endif
+                if ( maskoper.eq.'gt' ) then
+                   if ( mask(i,j).gt.maskvalue ) then
+                      nmask = nmask + 1
+                      mask(i,j) = 1.
+                   else
+                      mask(i,j) = 0.
+                   endif
+                endif
+                if ( maskoper.eq.'lt' ) then
+                   if ( mask(i,j).lt.maskvalue ) then
+                      nmask = nmask + 1
+                      mask(i,j) = 1.
+                   else
+                      mask(i,j) = 0.
+                   endif
+                endif
+                if ( maskoper.eq.'ge' ) then
+                   if ( mask(i,j).ge.maskvalue ) then
+                      nmask = nmask + 1
+                      mask(i,j) = 1.
+                   else
+                      mask(i,j) = 0.
+                   endif
+                endif
+                if ( maskoper.eq.'le' ) then
+                   if ( mask(i,j).le.maskvalue ) then
+                      nmask = nmask + 1
+                      mask(i,j) = 1.
+                   else
+                      mask(i,j) = 0.
+                   endif
+                endif
+             enddo
+          enddo
+          print*
+          print*,'  Mask read from file -> ',nmask,' 1-points'
+          print*
+      endif
 c     Write some status information
       print*,'---- CONSTANT GRID PARAMETERS ---------------------------'
       print*
       print*,'  xmin,xmax        : ',xmin,xmax
       print*,'  ymin,ymax        : ',ymin,ymax
                   lonlist(hn) = lon
                   latlist(hn) = lat
                endif
             enddo
          enddo
       endif
+c     mask: grid
+      if ( hmode.eq.'mask.grid' ) then
+         hn  = 0
+         do i=1,nx
+          do j=1,ny
+            if ( mask(i,j).gt.0.5 ) then
+               hn = hn +1
+               lonlist(hn) = xmin + real(i-1) * dx
+               latlist(hn) = ymin + real(j-1) * dy
+            endif
+          enddo
+         enddo
+      endif
+c     mask: equidistant
+      if ( hmode.eq.'mask.eqd' ) then
+         hn  = 0
+         lat = -90.+dy
+         do while ( lat.le.ymax )
+           lon = -180.
+           do while ( lon.le.xmax )
+              indx = nint( (lon-xmin)/dx + 1. )
+              indy = nint( (lat-ymin)/dy + 1. )
+              if ( indx.lt. 1 ) indx =  1
+              if ( indy.lt. 1 ) indy =  1
+              if ( indx.gt.nx ) indx = nx
+              if ( indy.gt.ny ) indy = ny
+              if ( mask(indx,indy).gt.0.5 ) then
+                 hn          = hn+1
+                 lonlist(hn) = lon
+                 latlist(hn) = lat
+              endif
+              lon = lon+ds/(deltat*cos(pi180*lat))
+           enddo
+           lat = lat+ds/deltat
+        enddo
+      endif
 c     ------ Get lev (vertical) coordinates -------------------------
 c     Read level list from file
       if ( vmode.eq.'file' ) then
          vn = 0

Subversion Repositories lagranto.ecmwf

(root)/trunk/create_startf/create_startf.f – Rev 15 → 19