WebSVN – lagranto.wrf – Diff – /trunk/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
      real,allocatable,dimension (:,:  )  :: temp3
      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 = 'nil'
          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) then
             print*,'Error in inq:',trim(varname)
         endif
         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) )
          allocate( temp3(nx,ny   ) )
 
      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 W"
 
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
                                            
c	     status = NF_INQ_VARID (fid, 'PHB', varid)
c	     if (status .ne. NF_NOERR) print*,"Error in inq sgeopot"
c	     status = NF_GET_VAR_REAL (fid, varid, temp1)
c	     if (status .ne. NF_NOERR) print*,"Error in reading sgeopot"
c
c	     status = NF_INQ_VARID (fid, 'PH', varid)
c	     if (status .ne. NF_NOERR) print*,"Error in inq pgeopot"
c	     status = NF_GET_VAR_REAL (fid, varid, temp2)
c	     if (status .ne. NF_NOERR) print*,"Error in reading pgeopot"
c
c	     do k = 1, nz+1
c	       do j = 1, ny
c	         do i = 1, nx
c	           temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)
c	         enddo
c	       enddo
c	     enddo
 
	     status = NF_INQ_VARID (fid, 'HGT', varid)
	     if (status .ne. NF_NOERR) print*,"Error in inq HGT"
	     status = NF_GET_VAR_REAL (fid, varid, temp3)
	     if (status .ne. NF_NOERR) print*,"Error in reading HGT"
	     do k = 1, nz
	       do j = 1, ny
	         do i = 1, nx
	           temp(i,j,k) = temp3(i,j) * gearth
	         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     ------------------------------------------------------------	

Rev 15	Rev 21
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	`real,allocatable,dimension (:,: ) :: temp3`
110	`integer ndims`	111	`integer ndims`
111	`integer is2d`	112	`integer is2d`
112	`integer levid,lonid,latid`	113	`integer levid,lonid,latid`
113	`integer nx,ny,nz`	114	`integer nx,ny,nz`
114	`integer status,i,j,k`	115	`integer status,i,j,k`
115	`character*2 stag`	116	`character*2 stag`
116		117
117	`c ------------ Get grid info ---------------------------------------`	118	`c ------------ Get grid info ---------------------------------------`
118	`status = NF_INQ_DIMID(fid, 'south_north', LATID)`	119	`status = NF_INQ_DIMID(fid, 'south_north', LATID)`
119	`status = NF_INQ_DIMID(fid, 'west_east', LONID)`	120	`status = NF_INQ_DIMID(fid, 'west_east', LONID)`
120	`status = NF_INQ_DIMID(fid, 'bottom_top', LEVID)`	121	`status = NF_INQ_DIMID(fid, 'bottom_top', LEVID)`
121	`if (status .ne. NF_NOERR) print*,"Error in reading grid atts"`	122	`if (status .ne. NF_NOERR) print*,"Error in reading grid atts"`
122		123
123	`status = NF_INQ_DIMLEN(fid, LONID, nx)`	124	`status = NF_INQ_DIMLEN(fid, LONID, nx)`
124	`status = NF_INQ_DIMLEN(fid, LATID, ny)`	125	`status = NF_INQ_DIMLEN(fid, LATID, ny)`
125	`status = NF_INQ_DIMLEN(fid, LEVID, nz)`	126	`status = NF_INQ_DIMLEN(fid, LEVID, nz)`
126	`if (status .ne. NF_NOERR) print*,"Error in reading grid size"`	127	`if (status .ne. NF_NOERR) print*,"Error in reading grid size"`
127		128
128	`c ------------ Allocate memory and set staggering mode -------------`	129	`c ------------ Allocate memory and set staggering mode -------------`
129	`if (trim(varname).eq.'U') then`	130	`if (trim(varname).eq.'U') then`
130	`stag = 'X'`	131	`stag = 'X'`
131	`is2d = 0`	132	`is2d = 0`
132		133
133	`elseif (trim(varname).eq.'V') then`	134	`elseif (trim(varname).eq.'V') then`
134	`stag = 'Y'`	135	`stag = 'Y'`
135	`is2d = 0`	136	`is2d = 0`
136		137
137	`elseif (trim(varname).eq.'W') then`	138	`elseif (trim(varname).eq.'W') then`
138	`stag = 'Z'`	139	`stag = 'Z'`
139	`is2d = 0`	140	`is2d = 0`
140		141
141	`elseif ( (trim(varname).eq.'geopot').or.`	142	`elseif ( (trim(varname).eq.'geopot').or.`
142	`> (trim(varname).eq.'Z' ) )`	143	`> (trim(varname).eq.'Z' ) )`
143	`>then`	144	`>then`
144	`stag = 'Z'`	145	`stag = 'Z'`
145	`is2d = 0`	146	`is2d = 0`
146		147
147	`elseif ( (trim(varname).eq.'geopots').or.`	148	`elseif ( (trim(varname).eq.'geopots').or.`
148	`> (trim(varname).eq.'ZB' ) )`	149	`> (trim(varname).eq.'ZB' ) )`
149	`>then`	150	`>then`
150	`stag = 'Z'`	151	`stag = 'nil'`
151	`is2d = 1`	152	`is2d = 1`
152		153
153	`elseif ( (trim(varname).eq.'pressure').or.`	154	`elseif ( (trim(varname).eq.'pressure').or.`
154	`> (trim(varname).eq.'P' ) )`	155	`> (trim(varname).eq.'P' ) )`
155	`>then`	156	`>then`
156	`stag = 'nil'`	157	`stag = 'nil'`
157	`is2d = 0`	158	`is2d = 0`
158		159
159	`else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then`	160	`else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then`
160	`stag = 'nil'`	161	`stag = 'nil'`
161	`is2d = 0`	162	`is2d = 0`
162		163
163	`else`	164	`else`
164	`status = NF_INQ_VARID (fid, trim(varname), varid)`	165	`status = NF_INQ_VARID (fid, trim(varname), varid)`
165	`status = NF_GET_ATT_TEXT(fid,varid,'stagger',stag)`	166	`status = NF_GET_ATT_TEXT(fid,varid,'stagger',stag)`
166	`if (status .ne. NF_NOERR) then`	167	`if (status .ne. NF_NOERR) then`
167	`print*,'Error in inq:',trim(varname)`	168	`print*,'Error in inq:',trim(varname)`
168	`endif`	169	`endif`
169	`status = NF_INQ_VARNDIMS (fid, varid, ndims)`	170	`status = NF_INQ_VARNDIMS (fid, varid, ndims)`
170	`is2d = 0`	171	`is2d = 0`
171	`if ( ndims.eq.3 ) is2d = 1`	172	`if ( ndims.eq.3 ) is2d = 1`
172		173
173	`endif`	174	`endif`
174		175
175	`if ( stag.eq.'X' ) then`	176	`if ( stag.eq.'X' ) then`
176	`allocate( temp (nx+1,ny,nz) )`	177	`allocate( temp (nx+1,ny,nz) )`
177	`allocate( temp1(nx+1,ny,nz) )`	178	`allocate( temp1(nx+1,ny,nz) )`
178	`allocate( temp2(nx+1,ny,nz) )`	179	`allocate( temp2(nx+1,ny,nz) )`
179		180
180	`elseif ( stag.eq.'Y' ) then`	181	`elseif ( stag.eq.'Y' ) then`
181	`allocate( temp (nx,ny+1,nz) )`	182	`allocate( temp (nx,ny+1,nz) )`
182	`allocate( temp1(nx,ny+1,nz) )`	183	`allocate( temp1(nx,ny+1,nz) )`
183	`allocate( temp2(nx,ny+1,nz) )`	184	`allocate( temp2(nx,ny+1,nz) )`
184		185
185	`elseif ( stag.eq.'Z' ) then`	186	`elseif ( stag.eq.'Z' ) then`
186	`allocate( temp (nx,ny,nz+1) )`	187	`allocate( temp (nx,ny,nz+1) )`
187	`allocate( temp1(nx,ny,nz+1) )`	188	`allocate( temp1(nx,ny,nz+1) )`
188	`allocate( temp2(nx,ny,nz+1) )`	189	`allocate( temp2(nx,ny,nz+1) )`
189		-
190	`else`	190	`else`
191	`allocate( temp (nx,ny,nz) )`	191	`allocate( temp (nx,ny,nz) )`
192	`allocate( temp1(nx,ny,nz) )`	192	`allocate( temp1(nx,ny,nz) )`
193	`allocate( temp2(nx,ny,nz) )`	193	`allocate( temp2(nx,ny,nz) )`
-		194	`allocate( temp3(nx,ny ) )`
194		195
195	`endif`	196	`endif`
196		197
197	`c ------------ Read data ------------------------------------------`	198	`c ------------ Read data ------------------------------------------`
198		199
199	`c Zonal wind : temp(nx+1,ny,nk)`	200	`c Zonal wind : temp(nx+1,ny,nk)`
200	`if (trim(varname).eq.'U') then`	201	`if (trim(varname).eq.'U') then`
201		202
202	`status = NF_INQ_VARID (fid, 'U', varid)`	203	`status = NF_INQ_VARID (fid, 'U', varid)`
203	`if (status .ne. NF_NOERR) print*,"Error in inq U"`	204	`if (status .ne. NF_NOERR) print*,"Error in inq U"`
204	`status = NF_GET_VAR_REAL (fid, varid, temp)`	205	`status = NF_GET_VAR_REAL (fid, varid, temp)`
205	`if (status .ne. NF_NOERR) print*,"Error in reading U"`	206	`if (status .ne. NF_NOERR) print*,"Error in reading U"`
206		207
207	`c Meridional wind : temp(nx,ny+1,nk)`	208	`c Meridional wind : temp(nx,ny+1,nk)`
208	`else if (trim(varname).eq.'V') then`	209	`else if (trim(varname).eq.'V') then`
209		210
210	`status = NF_INQ_VARID (fid, 'V', varid)`	211	`status = NF_INQ_VARID (fid, 'V', varid)`
211	`if (status .ne. NF_NOERR) print*,"Error in inq V"`	212	`if (status .ne. NF_NOERR) print*,"Error in inq V"`
212	`status = NF_GET_VAR_REAL (fid, varid, temp)`	213	`status = NF_GET_VAR_REAL (fid, varid, temp)`
213	`if (status .ne. NF_NOERR) print*,"Error in reading V"`	214	`if (status .ne. NF_NOERR) print*,"Error in reading V"`
214		215
215	`c Vertical wind : temp(nx,ny,nz+1)`	216	`c Vertical wind : temp(nx,ny,nz+1)`
216	`else if (trim(varname).eq.'W') then`	217	`else if (trim(varname).eq.'W') then`
217		218
218	`status = NF_INQ_VARID (fid, 'W', varid)`	219	`status = NF_INQ_VARID (fid, 'W', varid)`
219	`if (status .ne. NF_NOERR) print*,"Error in inq W"`	220	`if (status .ne. NF_NOERR) print*,"Error in inq W"`
220	`status = NF_GET_VAR_REAL (fid, varid, temp)`	221	`status = NF_GET_VAR_REAL (fid, varid, temp)`
221	`if (status .ne. NF_NOERR) print*,"Error in reading U"`	222	`if (status .ne. NF_NOERR) print*,"Error in reading W"`
222		223
223	`c Geopotential height (base + perturbation) : temp(nx,ny,nz+1)`	224	`c Geopotential height (base + perturbation) : temp(nx,ny,nz+1)`
224	`else if ( (trim(varname).eq.'geopot').or.`	225	`else if ( (trim(varname).eq.'geopot').or.`
225	`> (trim(varname).eq.'Z' ) )`	226	`> (trim(varname).eq.'Z' ) )`
226	`>then`	227	`>then`
227		228
228	`status = NF_INQ_VARID (fid, 'PHB', varid)`	229	`status = NF_INQ_VARID (fid, 'PHB', varid)`
229	`if (status .ne. NF_NOERR) print*,"Error in inq geopot"`	230	`if (status .ne. NF_NOERR) print*,"Error in inq geopot"`
230	`status = NF_GET_VAR_REAL (fid, varid, temp1)`	231	`status = NF_GET_VAR_REAL (fid, varid, temp1)`
231	`if (status .ne. NF_NOERR) print*,"Error in reading geopot"`	232	`if (status .ne. NF_NOERR) print*,"Error in reading geopot"`
232		233
233	`status = NF_INQ_VARID (fid, 'PH', varid)`	234	`status = NF_INQ_VARID (fid, 'PH', varid)`
234	`if (status .ne. NF_NOERR) print*,"Error in inq pgeopot"`	235	`if (status .ne. NF_NOERR) print*,"Error in inq pgeopot"`
235	`status = NF_GET_VAR_REAL (fid, varid, temp2)`	236	`status = NF_GET_VAR_REAL (fid, varid, temp2)`
236	`if (status .ne. NF_NOERR) print*,"Error in reading pgeopot"`	237	`if (status .ne. NF_NOERR) print*,"Error in reading pgeopot"`
237		238
238	`do k = 1, nz+1`	239	`do k = 1, nz+1`
239	`do j = 1, ny`	240	`do j = 1, ny`
240	`do i = 1, nx`	241	`do i = 1, nx`
241	`temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)`	242	`temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)`
242	`enddo`	243	`enddo`
243	`enddo`	244	`enddo`
244	`enddo`	245	`enddo`
245		246
246	`c surface geopotential: temp(nx,ny,nz+1)`	247	`c surface geopotential: temp(nx,ny,nz+1)`
247	`else if ( (trim(varname).eq.'geopots').or.`	248	`else if ( (trim(varname).eq.'geopots').or.`
248	`> (trim(varname).eq.'ZB' ) )`	249	`> (trim(varname).eq.'ZB' ) )`
249	`>then`	250	`>then`
250		251
251	`status = NF_INQ_VARID (fid, 'PHB', varid)`	252	`c status = NF_INQ_VARID (fid, 'PHB', varid)`
252	`if (status .ne. NF_NOERR) print*,"Error in inq sgeopot"`	253	`c if (status .ne. NF_NOERR) print*,"Error in inq sgeopot"`
253	`status = NF_GET_VAR_REAL (fid, varid, temp1)`	254	`c status = NF_GET_VAR_REAL (fid, varid, temp1)`
254	`if (status .ne. NF_NOERR) print*,"Error in reading sgeopot"`	255	`c if (status .ne. NF_NOERR) print*,"Error in reading sgeopot"`
255		256	`c`
256	`status = NF_INQ_VARID (fid, 'PH', varid)`	257	`c status = NF_INQ_VARID (fid, 'PH', varid)`
257	`if (status .ne. NF_NOERR) print*,"Error in inq pgeopot"`	258	`c if (status .ne. NF_NOERR) print*,"Error in inq pgeopot"`
258	`status = NF_GET_VAR_REAL (fid, varid, temp2)`	259	`c status = NF_GET_VAR_REAL (fid, varid, temp2)`
259	`if (status .ne. NF_NOERR) print*,"Error in reading pgeopot"`	260	`c if (status .ne. NF_NOERR) print*,"Error in reading pgeopot"`
-		261	`c`
-		262	`c do k = 1, nz+1`
-		263	`c do j = 1, ny`
-		264	`c do i = 1, nx`
-		265	`c temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)`
-		266	`c enddo`
-		267	`c enddo`
-		268	`c enddo`
260		269
-		270	`status = NF_INQ_VARID (fid, 'HGT', varid)`
-		271	`if (status .ne. NF_NOERR) print*,"Error in inq HGT"`
-		272	`status = NF_GET_VAR_REAL (fid, varid, temp3)`
-		273	`if (status .ne. NF_NOERR) print*,"Error in reading HGT"`
261	`do k = 1, nz+1`	274	`do k = 1, nz`
262	`do j = 1, ny`	275	`do j = 1, ny`
263	`do i = 1, nx`	276	`do i = 1, nx`
264	`temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)`	277	`temp(i,j,k) = temp3(i,j) * gearth`
265	`enddo`	278	`enddo`
266	`enddo`	279	`enddo`
267	`enddo`	280	`enddo`
268		281
269	`c Pressure (base + perturbation) : temp(nx,ny,nz)`	282	`c Pressure (base + perturbation) : temp(nx,ny,nz)`
270	`elseif ( (trim(varname).eq.'pressure').or.`	283	`elseif ( (trim(varname).eq.'pressure').or.`
271	`> (trim(varname).eq.'P' ) )`	284	`> (trim(varname).eq.'P' ) )`
272	`>then`	285	`>then`
273		286
274	`status = NF_INQ_VARID (fid, 'PB', varid)`	287	`status = NF_INQ_VARID (fid, 'PB', varid)`
275	`if (status .ne. NF_NOERR) print*,"Error in inq pb"`	288	`if (status .ne. NF_NOERR) print*,"Error in inq pb"`
276	`status = NF_GET_VAR_REAL (fid, varid, temp1)`	289	`status = NF_GET_VAR_REAL (fid, varid, temp1)`
277	`if (status .ne. NF_NOERR) print*,"Error in reading pb"`	290	`if (status .ne. NF_NOERR) print*,"Error in reading pb"`
278		291
279	`status = NF_INQ_VARID (fid, 'P', varid)`	292	`status = NF_INQ_VARID (fid, 'P', varid)`
280	`if (status .ne. NF_NOERR) print*,"Error in inq p"`	293	`if (status .ne. NF_NOERR) print*,"Error in inq p"`
281	`status = NF_GET_VAR_REAL (fid, varid, temp2)`	294	`status = NF_GET_VAR_REAL (fid, varid, temp2)`
282	`if (status .ne. NF_NOERR) print*,"Error in reading p"`	295	`if (status .ne. NF_NOERR) print*,"Error in reading p"`
283		296
284	`do k = 1, nz`	297	`do k = 1, nz`
285	`do j = 1, ny`	298	`do j = 1, ny`
286	`do i = 1, nx`	299	`do i = 1, nx`
287	`temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)`	300	`temp(i,j,k) = temp1(i,j,k) + temp2(i,j,k)`
288	`enddo`	301	`enddo`
289	`enddo`	302	`enddo`
290	`enddo`	303	`enddo`
291		304
292	`c Potential temperature: temp(nx,ny,nz)`	305	`c Potential temperature: temp(nx,ny,nz)`
293	`else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then`	306	`else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then`
294		307
295	`status = NF_INQ_VARID (fid, 'T', varid)`	308	`status = NF_INQ_VARID (fid, 'T', varid)`
296	`if (status .ne. NF_NOERR) print*,"Error in inq T"`	309	`if (status .ne. NF_NOERR) print*,"Error in inq T"`
297	`status = NF_GET_VAR_REAL (fid, varid, temp)`	310	`status = NF_GET_VAR_REAL (fid, varid, temp)`
298	`if (status .ne. NF_NOERR) print*,"Error in reading T"`	311	`if (status .ne. NF_NOERR) print*,"Error in reading T"`
299		312
300	`c Any other field (2D + 3D): temp(:,:,:), depending on staggering`	313	`c Any other field (2D + 3D): temp(:,:,:), depending on staggering`
301	`else`	314	`else`
302		315
303	`status = NF_INQ_VARID (fid, trim(varname), varid)`	316	`status = NF_INQ_VARID (fid, trim(varname), varid)`
304	`if (status .ne. NF_NOERR) then`	317	`if (status .ne. NF_NOERR) then`
305	`print*,"Error in inq:",trim(varname)`	318	`print*,"Error in inq:",trim(varname)`
306	`endif`	319	`endif`
307	`status = NF_GET_VAR_REAL (fid, varid, temp)`	320	`status = NF_GET_VAR_REAL (fid, varid, temp)`
308	`if (status .ne. NF_NOERR) then`	321	`if (status .ne. NF_NOERR) then`
309	`print*,"Error in reading:",trim(varname)`	322	`print*,"Error in reading:",trim(varname)`
310	`endif`	323	`endif`
311		324
312	`endif`	325	`endif`
313		326
314	`c ------------ Destaggering in X, Y and Z direction ---------------`	327	`c ------------ Destaggering in X, Y and Z direction ---------------`
315	`if (trim(stag).eq.'X') then`	328	`if (trim(stag).eq.'X') then`
316	`do k = 1, n3`	329	`do k = 1, n3`
317	`do j = 1, n2`	330	`do j = 1, n2`
318	`do i = 1, n1`	331	`do i = 1, n1`
319	`field(i,j,k) = 0.5*(temp(i,j,k) + temp(i+1,j,k))`	332	`field(i,j,k) = 0.5*(temp(i,j,k) + temp(i+1,j,k))`
320	`enddo`	333	`enddo`
321	`enddo`	334	`enddo`
322	`enddo`	335	`enddo`
323		336
324	`elseif (trim(stag).eq.'Y') then`	337	`elseif (trim(stag).eq.'Y') then`
325	`do k = 1, n3`	338	`do k = 1, n3`
326	`do j = 1, n2`	339	`do j = 1, n2`
327	`do i = 1, n1`	340	`do i = 1, n1`
328	`field(i,j,k) = 0.5*(temp(i,j,k) + temp(i,j+1,k))`	341	`field(i,j,k) = 0.5*(temp(i,j,k) + temp(i,j+1,k))`
329	`enddo`	342	`enddo`
330	`enddo`	343	`enddo`
331	`enddo`	344	`enddo`
332		345
333	`elseif (trim(stag).eq.'Z') then`	346	`elseif (trim(stag).eq.'Z') then`
334	`do k = 1, n3`	347	`do k = 1, n3`
335	`do j = 1, n2`	348	`do j = 1, n2`
336	`do i = 1, n1`	349	`do i = 1, n1`
337	`field(i,j,k) = 0.5*(temp(i,j,k) + temp(i,j,k+1))`	350	`field(i,j,k) = 0.5*(temp(i,j,k) + temp(i,j,k+1))`
338	`enddo`	351	`enddo`
339	`enddo`	352	`enddo`
340	`enddo`	353	`enddo`
341		354
342	`else`	355	`else`
343	`do k = 1, n3`	356	`do k = 1, n3`
344	`do j = 1, n2`	357	`do j = 1, n2`
345	`do i = 1, n1`	358	`do i = 1, n1`
346	`field(i,j,k) = temp(i,j,k)`	359	`field(i,j,k) = temp(i,j,k)`
347	`enddo`	360	`enddo`
348	`enddo`	361	`enddo`
349	`enddo`	362	`enddo`
350		363
351	`endif`	364	`endif`
352		365
353	`c ---------- Change units -----------------------------------------`	366	`c ---------- Change units -----------------------------------------`
354		367
355	`c Pressure Pa -> hPa`	368	`c Pressure Pa -> hPa`
356	`if (trim(varname).eq.'pressure' .or. trim(varname).eq.'P') then`	369	`if (trim(varname).eq.'pressure' .or. trim(varname).eq.'P') then`
357	`do k = 1, n3`	370	`do k = 1, n3`
358	`do j = 1, n2`	371	`do j = 1, n2`
359	`do i = 1, n1`	372	`do i = 1, n1`
360	`field(i,j,k)=0.01 * field(i,j,k)`	373	`field(i,j,k)=0.01 * field(i,j,k)`
361	`enddo`	374	`enddo`
362	`enddo`	375	`enddo`
363	`enddo`	376	`enddo`
364		377
365	`c Potential temperature + 300`	378	`c Potential temperature + 300`
366	`else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then`	379	`else if (trim(varname).eq.'TH' .or. trim(varname).eq.'T') then`
367	`do k = 1, n3`	380	`do k = 1, n3`
368	`do j = 1, n2`	381	`do j = 1, n2`
369	`do i = 1, n1`	382	`do i = 1, n1`
370	`field(i,j,k)=field(i,j,k) + addtheta`	383	`field(i,j,k)=field(i,j,k) + addtheta`
371	`enddo`	384	`enddo`
372	`enddo`	385	`enddo`
373	`enddo`	386	`enddo`
374		387
375	`c Geopotential -> Geopotential Height (3D + Surface)`	388	`c Geopotential -> Geopotential Height (3D + Surface)`
376	`else if ( (trim(varname).eq.'geopot' ).or.`	389	`else if ( (trim(varname).eq.'geopot' ).or.`
377	`> (trim(varname).eq.'Z' ).or.`	390	`> (trim(varname).eq.'Z' ).or.`
378	`> (trim(varname).eq.'geopots').or.`	391	`> (trim(varname).eq.'geopots').or.`
379	`> (trim(varname).eq.'ZB' ) )`	392	`> (trim(varname).eq.'ZB' ) )`
380	`>then`	393	`>then`
381	`do k = 1, n3`	394	`do k = 1, n3`
382	`do j = 1, n2`	395	`do j = 1, n2`
383	`do i = 1, n1`	396	`do i = 1, n1`
384	`field(i,j,k)=field(i,j,k)/gearth`	397	`field(i,j,k)=field(i,j,k)/gearth`
385	`enddo`	398	`enddo`
386	`enddo`	399	`enddo`
387	`enddo`	400	`enddo`
388		401
389	`endif`	402	`endif`
390		403
391	`c ---------- Copy lowest level for 2d tracing ---------------------`	404	`c ---------- Copy lowest level for 2d tracing ---------------------`
392	`if ( is2d.eq.1 ) then`	405	`if ( is2d.eq.1 ) then`
393	`do i=1,n1`	406	`do i=1,n1`
394	`do j=1,n2`	407	`do j=1,n2`
395	`do k=1,n3`	408	`do k=1,n3`
396	`field(i,j,k) = field(i,j,1)`	409	`field(i,j,k) = field(i,j,1)`
397	`enddo`	410	`enddo`
398	`enddo`	411	`enddo`
399	`enddo`	412	`enddo`
400	`endif`	413	`endif`
401		414
402		415
403	`end`	416	`end`
404		417
405	`c ------------------------------------------------------------`	418	`c ------------------------------------------------------------`
406	`c Close input file`	419	`c Close input file`
407	`c ------------------------------------------------------------`	420	`c ------------------------------------------------------------`
408		421
409	`subroutine input_close(fid)`	422	`subroutine input_close(fid)`
410		423
411	`c Close the input file with file identifier <fid>.`	424	`c Close the input file with file identifier <fid>.`
412		425
413	`implicit none`	426	`implicit none`
414		427
415	`include 'netcdf.inc'`	428	`include 'netcdf.inc'`
416		429
417	`c Declaration of subroutine parameters`	430	`c Declaration of subroutine parameters`
418	`integer fid`	431	`integer fid`
419		432
420	`c Auxiliary variables`	433	`c Auxiliary variables`
421	`integer ierr`	434	`integer ierr`
422		435
423	`c Close file`	436	`c Close file`
424	`call ncclos(fid,ierr)`	437	`call ncclos(fid,ierr)`
425		438
426	`end`	439	`end`
427		440
428	`c ------------------------------------------------------------`	-
429	`c ------------------------------------------------------------`	441	`c ------------------------------------------------------------`

Subversion Repositories lagranto.wrf

(root)/trunk/lib/libwrf.f – Rev 15 → 21