WebSVN – lagranto.wrf – Diff – /branches/met.nps.monterey/lib/libwrf.f


c     ************************************************************
c     * This package provides input routines to read if grid     *
c     * informations or field from WRF output			 *
c     *		                                                 *
c     * Author: Sebastian Schemm, ETH, 2013 			 *
c     * V.1.0 only support for ideal channel implemented         *
c     *      for input_grid_wrf				         *
c     *								 *
c     ************************************************************
 
c     ------------------------------------------------------------
c     Open input file
c     ------------------------------------------------------------
 
      subroutine input_open (fid,filename)
 
c     Open the input file with filename <filename> and return the
c     file identifier <fid> for further reference. 
 
      implicit none
 
      include 'netcdf.inc'
 
c     Declaration of subroutine parameters
      integer      fid              ! File identifier
      character*80 filename         ! Filename
 
c     Declaration of auxiliary variables
      integer      ierr
 
c     Open IVE netcdf file
      fid = ncopn(filename,NCWRITE,ierr)
      if (ierr.ne.0) goto 900
      
c     Exception handling
      return
      
 900  print*,'Cannot open input file ',trim(filename)
      stop
 
      end
 
c     ------------------------------------------------------------
c     Read information about the grid
c     ------------------------------------------------------------
	
	  subroutine input_grid_wrf(fid,
     >				  xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz)
 
  	  implicit none
 
      include 'netcdf.inc'
 
c     Declaration of subroutine parameters
      integer      fid                 ! File identifier
	  integer	   latid,lonid,levid
      real         xmin,xmax,ymin,ymax ! Domain size
      real         dx,dy               ! Horizontal resolution
      integer      nx,ny,nz            ! Grid dimensions
	  integer	   status
 
 
C	  get grid info
	  status = NF_INQ_DIMID(fid, 'south_north', LATID)
	  status = NF_INQ_DIMID(fid, 'west_east',   LONID)
	  status = NF_INQ_DIMID(fid, 'bottom_top',  LEVID)
	  if (status .ne. NF_NOERR) print*,"Error in reading grid atts"
 
	  status = NF_INQ_DIMLEN(fid, LONID, nx)
	  status = NF_INQ_DIMLEN(fid, LATID, ny)
	  status = NF_INQ_DIMLEN(fid, LEVID, nz)
	  if (status .ne. NF_NOERR) print*,"Error in reading grid size"
 
	  status = NF_GET_ATT_REAL (fid, NF_GLOBAL, 'DX', dx)
	  status = NF_GET_ATT_REAL (fid, NF_GLOBAL, 'DY', dy)
	  if (status .ne. NF_NOERR) print*,"Error in reading dx, dy"
 
c	  setup a pseudo grid
      xmin = 0.0
      ymin = 0.0
	  xmax = xmin+real(nx-1)*dx
      ymax = ymin+real(ny-1)*dy
 
	  end
 
c     ------------------------------------------------------------
c     Read variables
c     ------------------------------------------------------------
 
      subroutine input_var_wrf(fid,varname,field,n1,n2,n3)
 
      implicit none
 
      include 'netcdf.inc'
 
c     Declaration of subroutine parameters
      integer		  fid, varid
      integer      	  n1,n2,n3
      character*80	  varname
      real            field(n1,n2,n3)
 
c     Parmeters
      real            gearth         ! Earth's acceleration
      parameter       (gearth=9.81)
      real            addtheta       ! Offset for potential temperature
      parameter       (addtheta=300.)
 
c     Auxiliary variables
      real,allocatable,dimension (:,:,:)  :: temp,temp1,temp2
      integer         ndims
      integer         is2d
      integer		  levid,lonid,latid
      integer      	  nx,ny,nz
      integer		  status,i,j,k
      character*2	  stag
 
c     ------------ Get grid info ---------------------------------------
      status = NF_INQ_DIMID(fid, 'south_north', LATID)
      status = NF_INQ_DIMID(fid, 'west_east',   LONID)
      status = NF_INQ_DIMID(fid, 'bottom_top',  LEVID)
      if (status .ne. NF_NOERR) print*,"Error in reading grid atts"
 
      status = NF_INQ_DIMLEN(fid, LONID, nx)
      status = NF_INQ_DIMLEN(fid, LATID, ny)
      status = NF_INQ_DIMLEN(fid, LEVID, nz)
      if (status .ne. NF_NOERR) print*,"Error in reading grid size"
 
c     ------------ Allocate memory and set staggering mode -------------
      if (trim(varname).eq.'U') then
          stag = 'X'
          is2d = 0
 
      elseif (trim(varname).eq.'V') then
          stag = 'Y'
          is2d = 0
 
      elseif (trim(varname).eq.'W') then
          stag = 'Z'
          is2d = 0
 
      elseif ( (trim(varname).eq.'geopot').or.
     >         (trim(varname).eq.'Z'     ) )
     >then
          stag = 'Z'
          is2d = 0
 
      elseif ( (trim(varname).eq.'geopots').or.
     >         (trim(varname).eq.'ZB'     ) )
     >then
          stag = 'Z'
          is2d = 1
 
      elseif ( (trim(varname).eq.'pressure').or.
     >         (trim(varname).eq.'P'       ) )
     >then
	      stag = 'nil'
          is2d = 0
 
	  else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then
          stag = 'nil'
          is2d = 0
 
      else
         status = NF_INQ_VARID (fid, trim(varname), varid)
	     status = NF_GET_ATT_TEXT(fid,varid,'stagger',stag)
	     if (status .ne. NF_NOERR) print*,"Error in inq:",trim(varname)
         status = NF_INQ_VARNDIMS (fid, varid, ndims)
         is2d = 0
         if ( ndims.eq.3 ) is2d = 1
 
      endif
 
      if ( stag.eq.'X' ) then
          allocate( temp (nx+1,ny,nz) )
          allocate( temp1(nx+1,ny,nz) )
          allocate( temp2(nx+1,ny,nz) )
 
      elseif ( stag.eq.'Y' ) then
          allocate( temp (nx,ny+1,nz) )
          allocate( temp1(nx,ny+1,nz) )
          allocate( temp2(nx,ny+1,nz) )
 
      elseif ( stag.eq.'Z' ) then
          allocate( temp (nx,ny,nz+1) )
          allocate( temp1(nx,ny,nz+1) )
          allocate( temp2(nx,ny,nz+1) )
 
      else
          allocate( temp (nx,ny,nz) )
          allocate( temp1(nx,ny,nz) )
          allocate( temp2(nx,ny,nz) )
 
      endif
 
c	  ------------ Read data ------------------------------------------
 
c	  Zonal wind : temp(nx+1,ny,nk)
	  if (trim(varname).eq.'U') then
                     
	     status = NF_INQ_VARID (fid, 'U', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq U"
	     status = NF_GET_VAR_REAL (fid, varid, temp)
	     if (status .ne. NF_NOERR) print*,"Error in reading U"
 
c     Meridional wind : temp(nx,ny+1,nk)
	  else if (trim(varname).eq.'V') then
                                             
	     status = NF_INQ_VARID (fid, 'V', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq V"
	     status = NF_GET_VAR_REAL (fid, varid, temp)
	     if (status .ne. NF_NOERR) print*,"Error in reading V"
 
c	  Vertical wind : temp(nx,ny,nz+1)
	  else if (trim(varname).eq.'W') then
                                            
	     status = NF_INQ_VARID (fid, 'W', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq W"
	     status = NF_GET_VAR_REAL (fid, varid, temp)
	     if (status .ne. NF_NOERR) print*,"Error in reading U"
 
c	  Geopotential height (base + perturbation) : temp(nx,ny,nz+1)
	  else if ( (trim(varname).eq.'geopot').or.
     >          (trim(varname).eq.'Z'     ) )
     >then
                                              
	     status = NF_INQ_VARID (fid, 'PHB', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq geopot"
	     status = NF_GET_VAR_REAL (fid, varid, temp1)
	     if (status .ne. NF_NOERR) print*,"Error in reading geopot"
 
	     status = NF_INQ_VARID (fid, 'PH', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq pgeopot"
	     status = NF_GET_VAR_REAL (fid, varid, temp2)
	     if (status .ne. NF_NOERR) print*,"Error in reading pgeopot"
 
	     do k = 1, nz+1
	       do j = 1, ny
	         do i = 1, nx
	           temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)
	         enddo
	       enddo
	     enddo
 
c	  surface geopotential: temp(nx,ny,nz+1)
	  else if ( (trim(varname).eq.'geopots').or.
     >          (trim(varname).eq.'ZB'     ) )
     >then
                                            
	     status = NF_INQ_VARID (fid, 'PHB', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq sgeopot"
	     status = NF_GET_VAR_REAL (fid, varid, temp1)
	     if (status .ne. NF_NOERR) print*,"Error in reading sgeopot"
 
	     status = NF_INQ_VARID (fid, 'PH', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq pgeopot"
	     status = NF_GET_VAR_REAL (fid, varid, temp2)
	     if (status .ne. NF_NOERR) print*,"Error in reading pgeopot"
 
	     do k = 1, nz+1
	       do j = 1, ny
	         do i = 1, nx
	           temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)
	         enddo
	       enddo
	     enddo
 
c	  Pressure (base + perturbation) : temp(nx,ny,nz)
	  elseif ( (trim(varname).eq.'pressure').or.
     >         (trim(varname).eq.'P'       ) )
     >then
                               
	      status = NF_INQ_VARID (fid, 'PB', varid)
	      if (status .ne. NF_NOERR) print*,"Error in inq pb"
	      status = NF_GET_VAR_REAL (fid, varid, temp1)
	      if (status .ne. NF_NOERR) print*,"Error in reading pb"
 
       	  status = NF_INQ_VARID (fid, 'P', varid)
	      if (status .ne. NF_NOERR) print*,"Error in inq p"
	      status = NF_GET_VAR_REAL (fid, varid, temp2)
	      if (status .ne. NF_NOERR) print*,"Error in reading p"
 
	      do k = 1, nz
	        do j = 1, ny
	          do i = 1, nx
	            temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)
	          enddo
	        enddo
	      enddo
	
c     Potential temperature: temp(nx,ny,nz)
	  else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then
                    
	      status = NF_INQ_VARID (fid, 'T', varid)
	      if (status .ne. NF_NOERR) print*,"Error in inq T"
	      status = NF_GET_VAR_REAL (fid, varid, temp)
	      if (status .ne. NF_NOERR) print*,"Error in reading T"
 
c	  Any other field (2D + 3D): temp(:,:,:), depending on staggering
	  else
 
	      status = NF_INQ_VARID (fid, trim(varname), varid)
	      if (status .ne. NF_NOERR) then
	         print*,"Error in inq:",trim(varname)
	 	  endif
	 	  status = NF_GET_VAR_REAL (fid, varid, temp)
	 	  if (status .ne. NF_NOERR) then
	 	     print*,"Error in reading:",trim(varname)
          endif
 
      endif
 
c     ------------ Destaggering in X, Y and Z direction ---------------
      if (trim(stag).eq.'X') then
         do k = 1, n3
	  	 do j = 1, n2
	   	 do i = 1, n1
	    	 field(i,j,k) = 0.5*(temp(i,j,k) + temp(i+1,j,k))
	   	 enddo
	  	 enddo
	 	 enddo
 
	  elseif (trim(stag).eq.'Y') then
	 	 do k = 1, n3
	  	 do j = 1, n2
	   	 do i = 1, n1
	    	 field(i,j,k) = 0.5*(temp(i,j,k) + temp(i,j+1,k))
	   	 enddo
	  	 enddo
	 	 enddo
 
	  elseif (trim(stag).eq.'Z') then
	     do k = 1, n3
	  	 do j = 1, n2
	   	 do i = 1, n1
	    	 field(i,j,k) = 0.5*(temp(i,j,k) + temp(i,j,k+1))
	   	 enddo
	  	 enddo
	 	 enddo
 
	  else
	  	 do k = 1, n3
	  	 do j = 1, n2
	   	 do i = 1, n1
	    	 field(i,j,k) = temp(i,j,k)
	   	 enddo
	  	 enddo
	 	 enddo
 
	  endif
 
c     ---------- Change units -----------------------------------------
 
c	  Pressure Pa -> hPa
	  if (trim(varname).eq.'pressure' .or. trim(varname).eq.'P') then
	    do k = 1, n3
	      do j = 1, n2
	        do i = 1, n1
	          field(i,j,k)=0.01 * field(i,j,k)
	        enddo
	      enddo
	    enddo
 
c     Potential temperature + 300
      else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then
	    do k = 1, n3
	      do j = 1, n2
	        do i = 1, n1
	          field(i,j,k)=field(i,j,k) + addtheta
	        enddo
	      enddo
	    enddo
 
c     Geopotential -> Geopotential Height (3D + Surface)
	  else if ( (trim(varname).eq.'geopot' ).or.
     >          (trim(varname).eq.'Z'      ).or.
     >          (trim(varname).eq.'geopots').or.
     >          (trim(varname).eq.'ZB'     ) )
     >then
     	do k = 1, n3
	      do j = 1, n2
	        do i = 1, n1
	          field(i,j,k)=field(i,j,k)/gearth
	        enddo
	      enddo
	    enddo
 
	  endif
 
c     ---------- Copy lowest level for 2d tracing ---------------------
      if ( is2d.eq.1 ) then
         do i=1,n1
           do j=1,n2
             do k=1,n3
                field(i,j,k) = field(i,j,1)
             enddo
           enddo
         enddo
      endif
 
 
	  end
 
c     ------------------------------------------------------------
c     Close input file
c     ------------------------------------------------------------
 
      subroutine input_close(fid)
 
c     Close the input file with file identifier <fid>.
 
      implicit none
 
      include 'netcdf.inc'
 
c     Declaration of subroutine parameters
      integer fid
 
c     Auxiliary variables
      integer ierr
 
c     Close file
      call ncclos(fid,ierr)
 
      end
 
c     ------------------------------------------------------------
c     ------------------------------------------------------------	

Rev 2	Rev 10
1	`c ************************************************************`	1	`c ************************************************************`
2	`c * This package provides input routines to read if grid *`	2	`c * This package provides input routines to read if grid *`
3	`c * informations or field from WRF output *`	3	`c * informations or field from WRF output *`
4	`c * *`	4	`c * *`
5	`c * Author: Sebastian Schemm, ETH, 2013 *`	5	`c * Author: Sebastian Schemm, ETH, 2013 *`
6	`c * V.1.0 only support for ideal channel implemented *`	6	`c * V.1.0 only support for ideal channel implemented *`
7	`c * for input_grid_wrf *`	7	`c * for input_grid_wrf *`
8	`c * *`	8	`c * *`
9	`c ************************************************************`	9	`c ************************************************************`
10		10
11	`c ------------------------------------------------------------`	11	`c ------------------------------------------------------------`
12	`c Open input file`	12	`c Open input file`
13	`c ------------------------------------------------------------`	13	`c ------------------------------------------------------------`
14		14
15	`subroutine input_open (fid,filename)`	15	`subroutine input_open (fid,filename)`
16		16
17	`c Open the input file with filename <filename> and return the`	17	`c Open the input file with filename <filename> and return the`
18	`c file identifier <fid> for further reference.`	18	`c file identifier <fid> for further reference.`
19		19
20	`implicit none`	20	`implicit none`
21		21
22	`include 'netcdf.inc'`	22	`include 'netcdf.inc'`
23		23
24	`c Declaration of subroutine parameters`	24	`c Declaration of subroutine parameters`
25	`integer fid ! File identifier`	25	`integer fid ! File identifier`
26	`character*80 filename ! Filename`	26	`character*80 filename ! Filename`
27		27
28	`c Declaration of auxiliary variables`	28	`c Declaration of auxiliary variables`
29	`integer ierr`	29	`integer ierr`
30		30
31	`c Open IVE netcdf file`	31	`c Open IVE netcdf file`
32	`fid = ncopn(filename,NCWRITE,ierr)`	32	`fid = ncopn(filename,NCWRITE,ierr)`
33	`if (ierr.ne.0) goto 900`	33	`if (ierr.ne.0) goto 900`
34		34
35	`c Exception handling`	35	`c Exception handling`
36	`return`	36	`return`
37		37
38	`900 print*,'Cannot open input file ',trim(filename)`	38	`900 print*,'Cannot open input file ',trim(filename)`
39	`stop`	39	`stop`
40		40
41	`end`	41	`end`
42		42
43	`c ------------------------------------------------------------`	43	`c ------------------------------------------------------------`
44	`c Read information about the grid`	44	`c Read information about the grid`
45	`c ------------------------------------------------------------`	45	`c ------------------------------------------------------------`
46		46
47	`subroutine input_grid_wrf(fid,`	47	`subroutine input_grid_wrf(fid,`
48	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz)`	48	`> xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz)`
49		49
50	`implicit none`	50	`implicit none`
51		51
52	`include 'netcdf.inc'`	52	`include 'netcdf.inc'`
53		53
54	`c Declaration of subroutine parameters`	54	`c Declaration of subroutine parameters`
55	`integer fid ! File identifier`	55	`integer fid ! File identifier`
56	`integer latid,lonid,levid`	56	`integer latid,lonid,levid`
57	`real xmin,xmax,ymin,ymax ! Domain size`	57	`real xmin,xmax,ymin,ymax ! Domain size`
58	`real dx,dy ! Horizontal resolution`	58	`real dx,dy ! Horizontal resolution`
59	`integer nx,ny,nz ! Grid dimensions`	59	`integer nx,ny,nz ! Grid dimensions`
60	`integer status`	60	`integer status`
61		61
62		62
63	`C get grid info`	63	`C get grid info`
64	`status = NF_INQ_DIMID(fid, 'south_north', LATID)`	64	`status = NF_INQ_DIMID(fid, 'south_north', LATID)`
65	`status = NF_INQ_DIMID(fid, 'west_east', LONID)`	65	`status = NF_INQ_DIMID(fid, 'west_east', LONID)`
66	`status = NF_INQ_DIMID(fid, 'bottom_top', LEVID)`	66	`status = NF_INQ_DIMID(fid, 'bottom_top', LEVID)`
67	`if (status .ne. NF_NOERR) print*,"Error in reading grid atts"`	67	`if (status .ne. NF_NOERR) print*,"Error in reading grid atts"`
68		68
69	`status = NF_INQ_DIMLEN(fid, LONID, nx)`	69	`status = NF_INQ_DIMLEN(fid, LONID, nx)`
70	`status = NF_INQ_DIMLEN(fid, LATID, ny)`	70	`status = NF_INQ_DIMLEN(fid, LATID, ny)`
71	`status = NF_INQ_DIMLEN(fid, LEVID, nz)`	71	`status = NF_INQ_DIMLEN(fid, LEVID, nz)`
72	`if (status .ne. NF_NOERR) print*,"Error in reading grid size"`	72	`if (status .ne. NF_NOERR) print*,"Error in reading grid size"`
73		73
74	`status = NF_GET_ATT_REAL (fid, NF_GLOBAL, 'DX', dx)`	74	`status = NF_GET_ATT_REAL (fid, NF_GLOBAL, 'DX', dx)`
75	`status = NF_GET_ATT_REAL (fid, NF_GLOBAL, 'DY', dy)`	75	`status = NF_GET_ATT_REAL (fid, NF_GLOBAL, 'DY', dy)`
76	`if (status .ne. NF_NOERR) print*,"Error in reading dx, dy"`	76	`if (status .ne. NF_NOERR) print*,"Error in reading dx, dy"`
77		77
78	`c setup a pseudo grid`	78	`c setup a pseudo grid`
79	`xmin = 0.0`	79	`xmin = 0.0`
80	`ymin = 0.0`	80	`ymin = 0.0`
81	`xmax = xmin+real(nx-1)*dx`	81	`xmax = xmin+real(nx-1)*dx`
82	`ymax = ymin+real(ny-1)*dy`	82	`ymax = ymin+real(ny-1)*dy`
83		83
84	`end`	84	`end`
85		85
86	`c ------------------------------------------------------------`	86	`c ------------------------------------------------------------`
87	`c Read variables`	87	`c Read variables`
88	`c ------------------------------------------------------------`	88	`c ------------------------------------------------------------`
89		89
90	`subroutine input_var_wrf(fid,varname,field,n1,n2,n3)`	90	`subroutine input_var_wrf(fid,varname,field,n1,n2,n3)`
91		91
92	`implicit none`	92	`implicit none`
93		93
94	`include 'netcdf.inc'`	94	`include 'netcdf.inc'`
95		95
96	`c Declaration of subroutine parameters`	96	`c Declaration of subroutine parameters`
97	`integer fid, varid`	97	`integer fid, varid`
98	`integer n1,n2,n3`	98	`integer n1,n2,n3`
99	`character*80 varname`	99	`character*80 varname`
100	`real field(n1,n2,n3)`	100	`real field(n1,n2,n3)`
101		101
102	`c Parmeters`	102	`c Parmeters`
103	`real gearth ! Earth's acceleration`	103	`real gearth ! Earth's acceleration`
104	`parameter (gearth=9.81)`	104	`parameter (gearth=9.81)`
105	`real addtheta ! Offset for potential temperature`	105	`real addtheta ! Offset for potential temperature`
106	`parameter (addtheta=300.)`	106	`parameter (addtheta=300.)`
107		107
108	`c Auxiliary variables`	108	`c Auxiliary variables`
109	`real,allocatable,dimension (:,:,:) :: temp,temp1,temp2`	109	`real,allocatable,dimension (:,:,:) :: temp,temp1,temp2`
110	`integer ndims`	110	`integer ndims`
111	`integer is2d`	111	`integer is2d`
112	`integer levid,lonid,latid`	112	`integer levid,lonid,latid`
113	`integer nx,ny,nz`	113	`integer nx,ny,nz`
114	`integer status,i,j,k`	114	`integer status,i,j,k`
115	`character*2 stag`	115	`character*2 stag`
116		116
117	`c ------------ Get grid info ---------------------------------------`	117	`c ------------ Get grid info ---------------------------------------`
118	`status = NF_INQ_DIMID(fid, 'south_north', LATID)`	118	`status = NF_INQ_DIMID(fid, 'south_north', LATID)`
119	`status = NF_INQ_DIMID(fid, 'west_east', LONID)`	119	`status = NF_INQ_DIMID(fid, 'west_east', LONID)`
120	`status = NF_INQ_DIMID(fid, 'bottom_top', LEVID)`	120	`status = NF_INQ_DIMID(fid, 'bottom_top', LEVID)`
121	`if (status .ne. NF_NOERR) print*,"Error in reading grid atts"`	121	`if (status .ne. NF_NOERR) print*,"Error in reading grid atts"`
122		122
123	`status = NF_INQ_DIMLEN(fid, LONID, nx)`	123	`status = NF_INQ_DIMLEN(fid, LONID, nx)`
124	`status = NF_INQ_DIMLEN(fid, LATID, ny)`	124	`status = NF_INQ_DIMLEN(fid, LATID, ny)`
125	`status = NF_INQ_DIMLEN(fid, LEVID, nz)`	125	`status = NF_INQ_DIMLEN(fid, LEVID, nz)`
126	`if (status .ne. NF_NOERR) print*,"Error in reading grid size"`	126	`if (status .ne. NF_NOERR) print*,"Error in reading grid size"`
127		127
128	`c ------------ Allocate memory and set staggering mode -------------`	128	`c ------------ Allocate memory and set staggering mode -------------`
129	`if (trim(varname).eq.'U') then`	129	`if (trim(varname).eq.'U') then`
130	`stag = 'X'`	130	`stag = 'X'`
131	`is2d = 0`	131	`is2d = 0`
132		132
133	`elseif (trim(varname).eq.'V') then`	133	`elseif (trim(varname).eq.'V') then`
134	`stag = 'Y'`	134	`stag = 'Y'`
135	`is2d = 0`	135	`is2d = 0`
136		136
137	`elseif (trim(varname).eq.'W') then`	137	`elseif (trim(varname).eq.'W') then`
138	`stag = 'Z'`	138	`stag = 'Z'`
139	`is2d = 0`	139	`is2d = 0`
140		140
141	`elseif ( (trim(varname).eq.'geopot').or.`	141	`elseif ( (trim(varname).eq.'geopot').or.`
142	`> (trim(varname).eq.'Z' ) )`	142	`> (trim(varname).eq.'Z' ) )`
143	`>then`	143	`>then`
144	`stag = 'Z'`	144	`stag = 'Z'`
145	`is2d = 0`	145	`is2d = 0`
146		146
147	`elseif ( (trim(varname).eq.'geopots').or.`	147	`elseif ( (trim(varname).eq.'geopots').or.`
148	`> (trim(varname).eq.'ZB' ) )`	148	`> (trim(varname).eq.'ZB' ) )`
149	`>then`	149	`>then`
150	`stag = 'Z'`	150	`stag = 'Z'`
151	`is2d = 1`	151	`is2d = 1`
152		152
153	`elseif ( (trim(varname).eq.'pressure').or.`	153	`elseif ( (trim(varname).eq.'pressure').or.`
154	`> (trim(varname).eq.'P' ) )`	154	`> (trim(varname).eq.'P' ) )`
155	`>then`	155	`>then`
156	`stag = 'nil'`	156	`stag = 'nil'`
157	`is2d = 0`	157	`is2d = 0`
158		158
159	`else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then`	159	`else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then`
160	`stag = 'nil'`	160	`stag = 'nil'`
161	`is2d = 0`	161	`is2d = 0`
162		162
163	`else`	163	`else`
164	`status = NF_INQ_VARID (fid, trim(varname), varid)`	164	`status = NF_INQ_VARID (fid, trim(varname), varid)`
165	`status = NF_GET_ATT_TEXT(fid,varid,'stagger',stag)`	165	`status = NF_GET_ATT_TEXT(fid,varid,'stagger',stag)`
166	`if (status .ne. NF_NOERR) print*,"Error in inq:",trim(varname)`	166	`if (status .ne. NF_NOERR) print*,"Error in inq:",trim(varname)`
167	`status = NF_INQ_VARNDIMS (fid, varid, ndims)`	167	`status = NF_INQ_VARNDIMS (fid, varid, ndims)`
168	`is2d = 0`	168	`is2d = 0`
169	`if ( ndims.eq.3 ) is2d = 1`	169	`if ( ndims.eq.3 ) is2d = 1`
170		170
171	`endif`	171	`endif`
172		172
173	`if ( stag.eq.'X' ) then`	173	`if ( stag.eq.'X' ) then`
174	`allocate( temp (nx+1,ny,nz) )`	174	`allocate( temp (nx+1,ny,nz) )`
175	`allocate( temp1(nx+1,ny,nz) )`	175	`allocate( temp1(nx+1,ny,nz) )`
176	`allocate( temp2(nx+1,ny,nz) )`	176	`allocate( temp2(nx+1,ny,nz) )`
177		177
178	`elseif ( stag.eq.'Y' ) then`	178	`elseif ( stag.eq.'Y' ) then`
179	`allocate( temp (nx,ny+1,nz) )`	179	`allocate( temp (nx,ny+1,nz) )`
180	`allocate( temp1(nx,ny+1,nz) )`	180	`allocate( temp1(nx,ny+1,nz) )`
181	`allocate( temp2(nx,ny+1,nz) )`	181	`allocate( temp2(nx,ny+1,nz) )`
182		182
183	`elseif ( stag.eq.'Z' ) then`	183	`elseif ( stag.eq.'Z' ) then`
184	`allocate( temp (nx,ny,nz+1) )`	184	`allocate( temp (nx,ny,nz+1) )`
185	`allocate( temp1(nx,ny,nz+1) )`	185	`allocate( temp1(nx,ny,nz+1) )`
186	`allocate( temp2(nx,ny,nz+1) )`	186	`allocate( temp2(nx,ny,nz+1) )`
187		187
188	`else`	188	`else`
189	`allocate( temp (nx,ny,nz) )`	189	`allocate( temp (nx,ny,nz) )`
190	`allocate( temp1(nx,ny,nz) )`	190	`allocate( temp1(nx,ny,nz) )`
191	`allocate( temp2(nx,ny,nz) )`	191	`allocate( temp2(nx,ny,nz) )`
192		192
193	`endif`	193	`endif`
194		194
195	`c ------------ Read data ------------------------------------------`	195	`c ------------ Read data ------------------------------------------`
196		196
197	`c Zonal wind : temp(nx+1,ny,nk)`	197	`c Zonal wind : temp(nx+1,ny,nk)`
198	`if (trim(varname).eq.'U') then`	198	`if (trim(varname).eq.'U') then`
199		199
200	`status = NF_INQ_VARID (fid, 'U', varid)`	200	`status = NF_INQ_VARID (fid, 'U', varid)`
201	`if (status .ne. NF_NOERR) print*,"Error in inq U"`	201	`if (status .ne. NF_NOERR) print*,"Error in inq U"`
202	`status = NF_GET_VAR_REAL (fid, varid, temp)`	202	`status = NF_GET_VAR_REAL (fid, varid, temp)`
203	`if (status .ne. NF_NOERR) print*,"Error in reading U"`	203	`if (status .ne. NF_NOERR) print*,"Error in reading U"`
204		204
205	`c Meridional wind : temp(nx,ny+1,nk)`	205	`c Meridional wind : temp(nx,ny+1,nk)`
206	`else if (trim(varname).eq.'V') then`	206	`else if (trim(varname).eq.'V') then`
207		207
208	`status = NF_INQ_VARID (fid, 'V', varid)`	208	`status = NF_INQ_VARID (fid, 'V', varid)`
209	`if (status .ne. NF_NOERR) print*,"Error in inq V"`	209	`if (status .ne. NF_NOERR) print*,"Error in inq V"`
210	`status = NF_GET_VAR_REAL (fid, varid, temp)`	210	`status = NF_GET_VAR_REAL (fid, varid, temp)`
211	`if (status .ne. NF_NOERR) print*,"Error in reading V"`	211	`if (status .ne. NF_NOERR) print*,"Error in reading V"`
212		212
213	`c Vertical wind : temp(nx,ny,nz+1)`	213	`c Vertical wind : temp(nx,ny,nz+1)`
214	`else if (trim(varname).eq.'W') then`	214	`else if (trim(varname).eq.'W') then`
215		215
216	`status = NF_INQ_VARID (fid, 'W', varid)`	216	`status = NF_INQ_VARID (fid, 'W', varid)`
217	`if (status .ne. NF_NOERR) print*,"Error in inq W"`	217	`if (status .ne. NF_NOERR) print*,"Error in inq W"`
218	`status = NF_GET_VAR_REAL (fid, varid, temp)`	218	`status = NF_GET_VAR_REAL (fid, varid, temp)`
219	`if (status .ne. NF_NOERR) print*,"Error in reading U"`	219	`if (status .ne. NF_NOERR) print*,"Error in reading U"`
220		220
221	`c Geopotential height (base + perturbation) : temp(nx,ny,nz+1)`	221	`c Geopotential height (base + perturbation) : temp(nx,ny,nz+1)`
222	`else if ( (trim(varname).eq.'geopot').or.`	222	`else if ( (trim(varname).eq.'geopot').or.`
223	`> (trim(varname).eq.'Z' ) )`	223	`> (trim(varname).eq.'Z' ) )`
224	`>then`	224	`>then`
225		225
226	`status = NF_INQ_VARID (fid, 'PHB', varid)`	226	`status = NF_INQ_VARID (fid, 'PHB', varid)`
227	`if (status .ne. NF_NOERR) print*,"Error in inq geopot"`	227	`if (status .ne. NF_NOERR) print*,"Error in inq geopot"`
228	`status = NF_GET_VAR_REAL (fid, varid, temp1)`	228	`status = NF_GET_VAR_REAL (fid, varid, temp1)`
229	`if (status .ne. NF_NOERR) print*,"Error in reading geopot"`	229	`if (status .ne. NF_NOERR) print*,"Error in reading geopot"`
230		230
231	`status = NF_INQ_VARID (fid, 'PH', varid)`	231	`status = NF_INQ_VARID (fid, 'PH', varid)`
232	`if (status .ne. NF_NOERR) print*,"Error in inq pgeopot"`	232	`if (status .ne. NF_NOERR) print*,"Error in inq pgeopot"`
233	`status = NF_GET_VAR_REAL (fid, varid, temp2)`	233	`status = NF_GET_VAR_REAL (fid, varid, temp2)`
234	`if (status .ne. NF_NOERR) print*,"Error in reading pgeopot"`	234	`if (status .ne. NF_NOERR) print*,"Error in reading pgeopot"`
235		235
236	`do k = 1, nz+1`	236	`do k = 1, nz+1`
237	`do j = 1, ny`	237	`do j = 1, ny`
238	`do i = 1, nx`	238	`do i = 1, nx`
239	`temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)`	239	`temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)`
240	`enddo`	240	`enddo`
241	`enddo`	241	`enddo`
242	`enddo`	242	`enddo`
243		243
244	`c surface geopotential: temp(nx,ny,nz+1)`	244	`c surface geopotential: temp(nx,ny,nz+1)`
245	`else if ( (trim(varname).eq.'geopots').or.`	245	`else if ( (trim(varname).eq.'geopots').or.`
246	`> (trim(varname).eq.'ZB' ) )`	246	`> (trim(varname).eq.'ZB' ) )`
247	`>then`	247	`>then`
248		248
249	`status = NF_INQ_VARID (fid, 'PHB', varid)`	249	`status = NF_INQ_VARID (fid, 'PHB', varid)`
250	`if (status .ne. NF_NOERR) print*,"Error in inq sgeopot"`	250	`if (status .ne. NF_NOERR) print*,"Error in inq sgeopot"`
251	`status = NF_GET_VAR_REAL (fid, varid, temp1)`	251	`status = NF_GET_VAR_REAL (fid, varid, temp1)`
252	`if (status .ne. NF_NOERR) print*,"Error in reading sgeopot"`	252	`if (status .ne. NF_NOERR) print*,"Error in reading sgeopot"`
253		253
254	`status = NF_INQ_VARID (fid, 'PH', varid)`	254	`status = NF_INQ_VARID (fid, 'PH', varid)`
255	`if (status .ne. NF_NOERR) print*,"Error in inq pgeopot"`	255	`if (status .ne. NF_NOERR) print*,"Error in inq pgeopot"`
256	`status = NF_GET_VAR_REAL (fid, varid, temp2)`	256	`status = NF_GET_VAR_REAL (fid, varid, temp2)`
257	`if (status .ne. NF_NOERR) print*,"Error in reading pgeopot"`	257	`if (status .ne. NF_NOERR) print*,"Error in reading pgeopot"`
258		258
259	`do k = 1, nz+1`	259	`do k = 1, nz+1`
260	`do j = 1, ny`	260	`do j = 1, ny`
261	`do i = 1, nx`	261	`do i = 1, nx`
262	`temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)`	262	`temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)`
263	`enddo`	263	`enddo`
264	`enddo`	264	`enddo`
265	`enddo`	265	`enddo`
266		266
267	`c Pressure (base + perturbation) : temp(nx,ny,nz)`	267	`c Pressure (base + perturbation) : temp(nx,ny,nz)`
268	`elseif ( (trim(varname).eq.'pressure').or.`	268	`elseif ( (trim(varname).eq.'pressure').or.`
269	`> (trim(varname).eq.'P' ) )`	269	`> (trim(varname).eq.'P' ) )`
270	`>then`	270	`>then`
271		271
272	`status = NF_INQ_VARID (fid, 'PB', varid)`	272	`status = NF_INQ_VARID (fid, 'PB', varid)`
273	`if (status .ne. NF_NOERR) print*,"Error in inq pb"`	273	`if (status .ne. NF_NOERR) print*,"Error in inq pb"`
274	`status = NF_GET_VAR_REAL (fid, varid, temp1)`	274	`status = NF_GET_VAR_REAL (fid, varid, temp1)`
275	`if (status .ne. NF_NOERR) print*,"Error in reading pb"`	275	`if (status .ne. NF_NOERR) print*,"Error in reading pb"`
276		276
277	`status = NF_INQ_VARID (fid, 'P', varid)`	277	`status = NF_INQ_VARID (fid, 'P', varid)`
278	`if (status .ne. NF_NOERR) print*,"Error in inq p"`	278	`if (status .ne. NF_NOERR) print*,"Error in inq p"`
279	`status = NF_GET_VAR_REAL (fid, varid, temp2)`	279	`status = NF_GET_VAR_REAL (fid, varid, temp2)`
280	`if (status .ne. NF_NOERR) print*,"Error in reading p"`	280	`if (status .ne. NF_NOERR) print*,"Error in reading p"`
281		281
282	`do k = 1, nz`	282	`do k = 1, nz`
283	`do j = 1, ny`	283	`do j = 1, ny`
284	`do i = 1, nx`	284	`do i = 1, nx`
285	`temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)`	285	`temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)`
286	`enddo`	286	`enddo`
287	`enddo`	287	`enddo`
288	`enddo`	288	`enddo`
289		289
290	`c Potential temperature: temp(nx,ny,nz)`	290	`c Potential temperature: temp(nx,ny,nz)`
291	`else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then`	291	`else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then`
292		292
293	`status = NF_INQ_VARID (fid, 'T', varid)`	293	`status = NF_INQ_VARID (fid, 'T', varid)`
294	`if (status .ne. NF_NOERR) print*,"Error in inq T"`	294	`if (status .ne. NF_NOERR) print*,"Error in inq T"`
295	`status = NF_GET_VAR_REAL (fid, varid, temp)`	295	`status = NF_GET_VAR_REAL (fid, varid, temp)`
296	`if (status .ne. NF_NOERR) print*,"Error in reading T"`	296	`if (status .ne. NF_NOERR) print*,"Error in reading T"`
297		297
298	`c Any other field (2D + 3D): temp(:,:,:), depending on staggering`	298	`c Any other field (2D + 3D): temp(:,:,:), depending on staggering`
299	`else`	299	`else`
300		300
301	`status = NF_INQ_VARID (fid, trim(varname), varid)`	301	`status = NF_INQ_VARID (fid, trim(varname), varid)`
302	`if (status .ne. NF_NOERR) then`	302	`if (status .ne. NF_NOERR) then`
303	`print*,"Error in inq:",trim(varname)`	303	`print*,"Error in inq:",trim(varname)`
304	`endif`	304	`endif`
305	`status = NF_GET_VAR_REAL (fid, varid, temp)`	305	`status = NF_GET_VAR_REAL (fid, varid, temp)`
306	`if (status .ne. NF_NOERR) then`	306	`if (status .ne. NF_NOERR) then`
307	`print*,"Error in reading:",trim(varname)`	307	`print*,"Error in reading:",trim(varname)`
308	`endif`	308	`endif`
309		309
310	`endif`	310	`endif`
311		311
312	`c ------------ Destaggering in X, Y and Z direction ---------------`	312	`c ------------ Destaggering in X, Y and Z direction ---------------`
313	`if (trim(stag).eq.'X') then`	313	`if (trim(stag).eq.'X') then`
314	`do k = 1, n3`	314	`do k = 1, n3`
315	`do j = 1, n2`	315	`do j = 1, n2`
316	`do i = 1, n1`	316	`do i = 1, n1`
317	`field(i,j,k) = 0.5*(temp(i,j,k) + temp(i+1,j,k))`	317	`field(i,j,k) = 0.5*(temp(i,j,k) + temp(i+1,j,k))`
318	`enddo`	318	`enddo`
319	`enddo`	319	`enddo`
320	`enddo`	320	`enddo`
321		321
322	`elseif (trim(stag).eq.'Y') then`	322	`elseif (trim(stag).eq.'Y') then`
323	`do k = 1, n3`	323	`do k = 1, n3`
324	`do j = 1, n2`	324	`do j = 1, n2`
325	`do i = 1, n1`	325	`do i = 1, n1`
326	`field(i,j,k) = 0.5*(temp(i,j,k) + temp(i,j+1,k))`	326	`field(i,j,k) = 0.5*(temp(i,j,k) + temp(i,j+1,k))`
327	`enddo`	327	`enddo`
328	`enddo`	328	`enddo`
329	`enddo`	329	`enddo`
330		330
331	`elseif (trim(stag).eq.'Z') then`	331	`elseif (trim(stag).eq.'Z') then`
332	`do k = 1, n3`	332	`do k = 1, n3`
333	`do j = 1, n2`	333	`do j = 1, n2`
334	`do i = 1, n1`	334	`do i = 1, n1`
335	`field(i,j,k) = 0.5*(temp(i,j,k) + temp(i,j,k+1))`	335	`field(i,j,k) = 0.5*(temp(i,j,k) + temp(i,j,k+1))`
336	`enddo`	336	`enddo`
337	`enddo`	337	`enddo`
338	`enddo`	338	`enddo`
339		339
340	`else`	340	`else`
341	`do k = 1, n3`	341	`do k = 1, n3`
342	`do j = 1, n2`	342	`do j = 1, n2`
343	`do i = 1, n1`	343	`do i = 1, n1`
344	`field(i,j,k) = temp(i,j,k)`	344	`field(i,j,k) = temp(i,j,k)`
345	`enddo`	345	`enddo`
346	`enddo`	346	`enddo`
347	`enddo`	347	`enddo`
348		348
349	`endif`	349	`endif`
350		350
351	`c ---------- Change units -----------------------------------------`	351	`c ---------- Change units -----------------------------------------`
352		352
353	`c Pressure Pa -> hPa`	353	`c Pressure Pa -> hPa`
354	`if (trim(varname).eq.'pressure' .or. trim(varname).eq.'P') then`	354	`if (trim(varname).eq.'pressure' .or. trim(varname).eq.'P') then`
355	`do k = 1, n3`	355	`do k = 1, n3`
356	`do j = 1, n2`	356	`do j = 1, n2`
357	`do i = 1, n1`	357	`do i = 1, n1`
358	`field(i,j,k)=0.01 * field(i,j,k)`	358	`field(i,j,k)=0.01 * field(i,j,k)`
359	`enddo`	359	`enddo`
360	`enddo`	360	`enddo`
361	`enddo`	361	`enddo`
362		362
363	`c Potential temperature + 300`	363	`c Potential temperature + 300`
364	`else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then`	364	`else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then`
365	`do k = 1, n3`	365	`do k = 1, n3`
366	`do j = 1, n2`	366	`do j = 1, n2`
367	`do i = 1, n1`	367	`do i = 1, n1`
368	`field(i,j,k)=field(i,j,k) + addtheta`	368	`field(i,j,k)=field(i,j,k) + addtheta`
369	`enddo`	369	`enddo`
370	`enddo`	370	`enddo`
371	`enddo`	371	`enddo`
372		372
373	`c Geopotential -> Geopotential Height (3D + Surface)`	373	`c Geopotential -> Geopotential Height (3D + Surface)`
374	`else if ( (trim(varname).eq.'geopot' ).or.`	374	`else if ( (trim(varname).eq.'geopot' ).or.`
375	`> (trim(varname).eq.'Z' ).or.`	375	`> (trim(varname).eq.'Z' ).or.`
376	`> (trim(varname).eq.'geopots').or.`	376	`> (trim(varname).eq.'geopots').or.`
377	`> (trim(varname).eq.'ZB' ) )`	377	`> (trim(varname).eq.'ZB' ) )`
378	`>then`	378	`>then`
379	`do k = 1, n3`	379	`do k = 1, n3`
380	`do j = 1, n2`	380	`do j = 1, n2`
381	`do i = 1, n1`	381	`do i = 1, n1`
382	`field(i,j,k)=field(i,j,k)/gearth`	382	`field(i,j,k)=field(i,j,k)/gearth`
383	`enddo`	383	`enddo`
384	`enddo`	384	`enddo`
385	`enddo`	385	`enddo`
386		386
387	`endif`	387	`endif`
388		388
389	`c ---------- Copy lowest level for 2d tracing ---------------------`	389	`c ---------- Copy lowest level for 2d tracing ---------------------`
390	`if ( is2d.eq.1 ) then`	390	`if ( is2d.eq.1 ) then`
391	`do i=1,n1`	391	`do i=1,n1`
392	`do j=1,n2`	392	`do j=1,n2`
393	`do k=1,n3`	393	`do k=1,n3`
394	`field(i,j,k) = field(i,j,1)`	394	`field(i,j,k) = field(i,j,1)`
395	`enddo`	395	`enddo`
396	`enddo`	396	`enddo`
397	`enddo`	397	`enddo`
398	`endif`	398	`endif`
399		399
400		400
401	`end`	401	`end`
402		402
403	`c ------------------------------------------------------------`	403	`c ------------------------------------------------------------`
404	`c Close input file`	404	`c Close input file`
405	`c ------------------------------------------------------------`	405	`c ------------------------------------------------------------`
406		406
407	`subroutine input_close(fid)`	407	`subroutine input_close(fid)`
408		408
409	`c Close the input file with file identifier <fid>.`	409	`c Close the input file with file identifier <fid>.`
410		410
411	`implicit none`	411	`implicit none`
412		412
413	`include 'netcdf.inc'`	413	`include 'netcdf.inc'`
414		414
415	`c Declaration of subroutine parameters`	415	`c Declaration of subroutine parameters`
416	`integer fid`	416	`integer fid`
417		417
418	`c Auxiliary variables`	418	`c Auxiliary variables`
419	`integer ierr`	419	`integer ierr`
420		420
421	`c Close file`	421	`c Close file`
422	`call ncclos(fid,ierr)`	422	`call ncclos(fid,ierr)`
423		423
424	`end`	424	`end`
425		425
426	`c ------------------------------------------------------------`	426	`c ------------------------------------------------------------`
427	`c ------------------------------------------------------------`	427	`c ------------------------------------------------------------`

Subversion Repositories lagranto.wrf

(root)/branches/met.nps.monterey/lib/libwrf.f – Rev 2 → 10