Subversion Repositories lagranto.ecmwf

Rev

Rev 13 | Go to most recent revision | Show entire file | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 13 Rev 21
Line 1... Line 1...
1 
      program getmima
 1 
      PROGRAM getmima
2 
C     ===============
 - 
 
3 
 
 2 
 
- 
 
 3 
c     ***********************************************************************
4 
C     Get the minimum and maximum of a variable either
 4 
c     * Get minimum and maximum value of a varaiable on pressure            *
5 
C     in the whole 3d data domain, or an a specified pressure or
 5 
c     * or isentropic surface                                               *
6 
C     theta level.
 - 
 
7 
C     The program is invoked by the shell-script getmima.
 6 
c     * Michael Sprenger / Spring, summer 2016                              *
8 
C
 - 
 
9 
C     August 96		H. Wernli
 7 
c     ***********************************************************************
10 
 
 - 
 
11 
      implicit none
 - 
 
12 
 
 8 
 
13 
      integer	nzmax,ntmax
 9 
      use netcdf
14 
      parameter(nzmax=100,ntmax=200)
 - 
 
15 
 
 10 
 
16 
      real,allocatable, dimension (:) :: sp,varf,field,varl,tt
 - 
 
17 
      integer   stat
 11 
      implicit none
18 
 
 12 
 
- 
 
 13 
c     ----------------------------------------------------------------------
19 
      real	time(ntmax),level
 14 
c     Declaration of parameters and variables
- 
 
 15 
c     ----------------------------------------------------------------------
- 
 
 16
20 
      character*80 cdfnam,cstnam,varnam,tvar,clev
 17 
c     Interpolation method
- 
 
 18 
      integer   ipom
21 
      character*20 vnam(100)
 19 
      parameter (ipom=0)
- 
 
 20
- 
 
 21 
c     Flag for timecheck
22 
      character*1  mode
 22 
      character*80 timecheck
23 
      integer	cdfid,cstid,ierr,ndim,vardim(4)
 23 
      parameter    (timecheck = 'no' )
- 
 
 24
- 
 
 25 
c     netCDF fields
24 
      real	dx,dy,mdv,varmin(4),varmax(4),stag(4)
 26 
      real,allocatable, dimension (:,:)   :: sp,varf
25 
      real      aklev(nzmax),bklev(nzmax),aklay(nzmax),bklay(nzmax),
 27 
      real,allocatable, dimension (:,:,:) :: field,varl,tt,p3
- 
 
 28 
      real,allocatable, dimension (:)     :: ak,bk
- 
 
 29 
      integer                                stat
- 
 
 30 
      character*80                           cdfnam,varnam
- 
 
 31 
      character*80                           vnam(200)
- 
 
 32 
      integer                                nvars
- 
 
 33 
      integer	                             cdfid,ierr
26 
     >          ak(nzmax),bk(nzmax)
 34 
      real	                                 mdv
- 
 
 35 
      real                                   stagz
- 
 
 36
27 
      logical	prelev
 37 
c     Grid description
- 
 
 38 
      real      xmin,xmax       ! Zonal grid extension
- 
 
 39 
      real      ymin,ymax       ! Meridional grid extension
28 
      integer	nx,ny,nz,ntimes,ipom,nvars
 40 
      integer   nx,ny,nz        ! Grid dimensions
- 
 
 41 
      real      dx,dy           ! Horizontal grid resolution
29 
      real	pollon,pollat,xphys,yphys
 42 
      real	    pollon,pollat   ! Pole position
- 
 
 43 
 
30 
      integer	i,n
 44 
c     Output variables
31 
      integer	imin,jmin,imax,jmax,minind,maxind,kmin,kmax,hmin,hmax
 45 
      real      min,max                      ! Minimum & maximum value
32 
      real	min,max,lonmin,latmin,lonmax,latmax,pmin,pmax,tmin,tmax
 46 
      real      lonmin,lonmax,latmin,latmax  ! Position of min & max
33 
 
 47 
 
- 
 
 48 
c     Auxiliary variables
- 
 
 49 
      integer        iargc
- 
 
 50 
      character*(80) arg
- 
 
 51 
      real           rd
- 
 
 52 
      integer        i,j
- 
 
 53 
      integer        minindx,minindy,maxindx,maxindy
- 
 
 54 
      character*80   tvar
- 
 
 55 
      character*1    mode
- 
 
 56 
      character*80   clev
- 
 
 57 
      real	         xphys,yphys
- 
 
 58 
      real	         level
- 
 
 59 
      integer        flag
- 
 
 60
- 
 
 61 
c     Externals
34 
      real      lmstolm,phstoph
 62 
      real      lmstolm,phstoph
- 
 
 63 
      external  lmstolm,phstoph
35 
 
 64 
 
36 
      integer   iargc
 65 
c     ----------------------------------------------------------------------
37 
      character*(80) arg
 66 
c     Preparation - argument handling, grid paraemters, memor allocation
38 
      integer   flag
 67 
c     ----------------------------------------------------------------------
39
68 
 
40 
c     check for sufficient requested arguments
 69 
c     Check for sufficient requested arguments
41 
      if (iargc().lt.2) then
 70 
      if (iargc().lt.2) then
42 
        print*,'USAGE: getmima filename var ',
 71 
        print*,'USAGE: getmima filename var ',
43 
     >         '[lev in the form Pval or Tval]'
 72 
     >         '[lev in the form Pval or Tval]'
44 
        call exit(1)
 73 
        call exit(1)
45 
      endif
 74 
      endif
46
75
47 
c     read and transform input
 76 
c     Read and transform input
48 
      call getarg(1,arg)
 77 
      call getarg(1,arg)
49 
      cdfnam=trim(arg)
 78 
      cdfnam=trim(arg)
50
- 
 
51 
      call getarg(2,arg)
 79 
      call getarg(2,arg)
52 
      varnam=trim(arg)
 80 
      varnam=trim(arg)
53 
 
 81 
 
54 
      if (iargc().eq.3) then
 82 
      if (iargc().eq.3) then
55 
        call getarg(3,arg)
 83 
        call getarg(3,arg)
Line 61... Line 89...
61 
      else
 89 
      else
62 
        mode='X'
 90 
        mode='X'
63 
        level=0.
 91 
        level=0.
64 
      endif
 92 
      endif
65 
 
 93 
 
66 
      prelev=.true.
 94 
c     Init level type and minimum,maximum
67
- 
 
68 
      min= 1.e19
 95 
      min= 1.e19
69 
      max=-1.e19
 96 
      max=-1.e19
70
97 
 
- 
 
 98 
c     Open netCDF file
71 
C     Open files and get infos about data domain
 99 
      call input_open (cdfid,cdfnam)
72
100 
 
73 
      call cdfopn(trim(cdfnam),cdfid,ierr)
 101 
c     Get list of variables on file
74 
      call getcfn(cdfid,cstnam,ierr)
 102 
      call input_getvars (cdfid,vnam,nvars)
75 
      call cdfopn(trim(cstnam),cstid,ierr)
 - 
 
76 
      call getdef(cdfid,trim(varnam),ndim,mdv,vardim,
 - 
 
77 
     >            varmin,varmax,stag,ierr)
 - 
 
78 
      if (ierr.ne.0) goto 920
 - 
 
79
103 
 
80 
C     Get the data array and the levels
 104 
C     Get infos about data domain - in particular: dimensions
81
- 
 
82 
      nx=vardim(1)
 105 
      nx       = 1
83 
      ny=vardim(2)
 106 
      ny       = 1
84 
      nz=vardim(3)
 107 
      nz       = 1
85 
      call gettimes(cdfid,time,ntimes,ierr)
 - 
 
86 
 
 - 
 
87 
      call getgrid(cstid,dx,dy,ierr)
 - 
 
88 
      call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
 108 
      call input_grid (-cdfid,varnam,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
89 
      call getpole(cstid,pollon,pollat,ierr)
 109 
     >               0.,pollon,pollat,rd,rd,nz,rd,rd,rd,stagz,timecheck)
90 
 
 110 
 
91 
C     Get memory for dynamic arrays
 111 
C     Get memory for dynamic arrays
92
- 
 
93 
      allocate(sp(nx*ny),stat=stat)
 112 
      allocate(sp(nx,ny),stat=stat)
94 
      if (stat.ne.0) print*,'*** error allocating array sp ***'
 113 
      if (stat.ne.0) print*,'*** error allocating array sp ***'
95 
      allocate(varf(nx*ny),stat=stat)
 114 
      allocate(varf(nx,ny),stat=stat)
96 
      if (stat.ne.0) print*,'*** error allocating array varf ***'
 115 
      if (stat.ne.0) print*,'*** error allocating array varf ***'
97 
      allocate(varl(nx*ny*nz),stat=stat)
 116 
      allocate(varl(nx,ny,nz),stat=stat)
98 
      if (stat.ne.0) print*,'*** error allocating array varl ***'
 117 
      if (stat.ne.0) print*,'*** error allocating array varl ***'
99 
      allocate(tt(nx*ny*nz),stat=stat)
 118 
      allocate(tt(nx,ny,nz),stat=stat)
100 
      if (stat.ne.0) print*,'*** error allocating array tt ***'
 119 
      if (stat.ne.0) print*,'*** error allocating array tt ***'
101 
      allocate(field(nx*ny*nz),stat=stat)
 120 
      allocate(field(nx,ny,nz),stat=stat)
102 
      if (stat.ne.0) print*,'*** error allocating array field ***'
 121 
      if (stat.ne.0) print*,'*** error allocating array field ***'
- 
 
 122 
      allocate(ak(nz),stat=stat)
- 
 
 123 
      if (stat.ne.0) print*,'*** error allocating array ak ***'
- 
 
 124 
      allocate(bk(nz),stat=stat)
- 
 
 125 
      if (stat.ne.0) print*,'*** error allocating array bk ***'
103
126
- 
 
 127 
c     ----------------------------------------------------------------------
104 
      do n=1,ntimes
 128 
c     Load data
- 
 
 129 
c     ----------------------------------------------------------------------
105
130 
 
- 
 
 131 
c     Get grid info for variable - non-constant part  (P, PS, AK, BK)
- 
 
 132 
      call input_grid
106 
      call getdat(cdfid,trim(varnam),time(n),0,field,ierr)
 133 
     >       (cdfid,varnam,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
- 
 
 134 
     >        0.,pollon,pollat,varl,sp,nz,ak,bk,stagz,timecheck)
107
135
- 
 
 136 
c     Load Variable
108 
C     Define the appropriate ak and bk-arrays
 137 
      call input_wind (cdfid,varnam,field,0.,stagz,mdv,
- 
 
 138 
     >                 xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
109
139
- 
 
 140 
c     In case of isentropic level, read temperature
110 
      if (stag(3).eq.0.) then
 141 
      if (mode.eq.'T') then
111 
        do i=1,nz
 142 
            tvar="T"
112 
          ak(i)=aklev(i)
 143 
            do i=1,nvars
- 
 
 144 
              if (vnam(i).eq."THETA") tvar="THETA"
113 
          bk(i)=bklev(i)
 145 
              if (vnam(i).eq."TH")    tvar="TH"
114 
        enddo
 146 
            enddo
115 
      else
 147 
            if (tvar.eq.'T') then
- 
 
 148 
                 call input_wind (cdfid,tvar,tt,0.,stagz,mdv,xmin,xmax,
- 
 
 149 
     >                            ymin,ymax,dx,dy,nx,ny,nz,timecheck)
- 
 
 150 
                 call pottemp(varl,tt,sp,nx,ny,nz,ak,bk)
116 
        do i=1,nz
 151 
            else
117 
          ak(i)=aklay(i)
 152 
                 call input_wind (cdfid,tvar,tt,0.,stagz,mdv,xmin,xmax,
118 
          bk(i)=bklay(i)
 153 
     >                            ymin,ymax,dx,dy,nx,ny,nz,timecheck)
119 
        enddo
 154 
            endif
120 
      endif
 155 
      endif
121 
 
 156 
 
122 
      do i=1,nz
 - 
 
123 
        if (bk(i).ne.0.) prelev=.false.
 - 
 
124 
      enddo
 - 
 
125 
 
 - 
 
126 
      if (.not.prelev) call getdat(cdfid,'PS',time(n),0,sp,ierr)
 157 
c     ----------------------------------------------------------------------
127 
 
 - 
 
128 
C     Determine name of "temperature variable"
 158 
c     Do interpolation and get minimum, maximum
129 
 
 - 
 
130 
      call getvars(cdfid,nvars,vnam,ierr)
 - 
 
131 
      tvar="T"
 - 
 
132 
      do i=1,nvars
 - 
 
133 
        if (vnam(i).eq."THETA") tvar="THETA"
 159 
c     ----------------------------------------------------------------------
134 
        if (vnam(i).eq."TH")    tvar="TH"
 - 
 
135 
      enddo
 - 
 
136 
 
 160 
 
137 
C     If required do interpolation on pressure or theta level
 161 
C     If required do interpolation on pressure or theta level
138
- 
 
139 
      if (mode.eq.'P') then
 162 
      if (mode.eq.'P') then
140 
        call pres(varl,sp,nx,ny,nz,ak,bk)
 - 
 
141 
        call vipo(field,varl,level,varf,nx,ny,nz,mdv,ipom)
 163 
        call vipo(field,varl,level,varf,nx,ny,nz,mdv,ipom)
142 
      else if (mode.eq.'T') then
 164 
      elseif (mode.eq.'T') then
- 
 
 165 
        call thipo(field,varl,level,varf,nx,ny,nz,mdv,ipom)
143 
        if (tvar.eq.'T') then
 166 
      elseif (mode.eq.'L') then
144 
          call getdat(cdfid,'T',time(n),0,tt,ierr)
 167 
        do i=1,nx
145 
          call pottemp(varl,tt,sp,nx,ny,nz,ak,bk)
 168 
           do j=1,ny
146 
          call thipo(field,varl,level,varf,nx,ny,nz,mdv,ipom)
 169 
               varf(i,j) = field(i,j,nint(level))
- 
 
 170 
           enddo
147 
        else
 171 
        enddo
148 
          call getdat(cdfid,tvar,time(n),0,varl,ierr)
 172 
      elseif (mode.eq.'X') then
- 
 
 173 
         do i=1,nx
- 
 
 174 
            do j=1,ny
149 
          call thipo(field,varl,level,varf,nx,ny,nz,mdv,ipom)
 175 
                varf(i,j) = field(i,j,1)
- 
 
 176 
            enddo
150 
        endif
 177 
         enddo
151 
      endif
 178 
      endif
152
179
153 
C     Determine minimum and maximum and calculate their location
 180 
C     Determine minimum and maximum
154
- 
 
155 
      if (mode.ne.'X') then
 181 
      do i=1,nx
156 
        do i=1,nx*ny
 182 
        do j=1,ny
157 
          if ((varf(i).ne.mdv).and.(varf(i).lt.min)) then
 183 
            if ((varf(i,j).ne.mdv).and.(varf(i,j).lt.min)) then
158 
            min=varf(i)
 184 
               min     = varf(i,j)
159 
            minind=i
 185 
               minindx = i
160 
            tmin=time(n)
 186 
               minindy = j
161 
          else if ((varf(i).ne.mdv).and.(varf(i).gt.max)) then
 187 
            elseif ((varf(i,j).ne.mdv).and.(varf(i,j).gt.max)) then
162 
            max=varf(i)
 188 
               max     = varf(i,j)
163 
            maxind=i
 189 
               maxindx = i
164 
            tmax=time(n)
 - 
 
165 
          endif
 - 
 
166 
        enddo
 - 
 
167 
        imin=mod(minind-1,nx)+1
 - 
 
168 
        jmin=int((minind-1)/nx)+1
 - 
 
169 
        imax=mod(maxind-1,nx)+1
 - 
 
170 
        jmax=int((maxind-1)/nx)+1
 - 
 
171 
        lonmin=varmin(1)+(imin-1)*dx
 - 
 
172 
        latmin=varmin(2)+(jmin-1)*dy
 - 
 
173 
        pmin=level
 - 
 
174 
        lonmax=varmin(1)+(imax-1)*dx
 - 
 
175 
        latmax=varmin(2)+(jmax-1)*dy
 - 
 
176 
        pmax=level
 - 
 
177 
      else
 - 
 
178 
        do i=1,nx*ny*nz
 - 
 
179 
          if ((field(i).ne.mdv).and.(field(i).lt.min)) then
 - 
 
180 
            min=field(i)
 - 
 
181 
            minind=i
 - 
 
182 
            tmin=time(n)
 - 
 
183 
          else if ((field(i).ne.mdv).and.(field(i).gt.max)) then
 - 
 
184 
            max=field(i)
 - 
 
185 
            maxind=i
 190 
               maxindy = j
186 
            tmax=time(n)
 - 
 
187 
          endif
 191 
          endif
188 
        enddo
 192 
        enddo
189 
        imin=mod(mod(minind-1,nx*ny),nx)+1
 - 
 
190 
        jmin=int(mod(minind-1,nx*ny)/nx)+1
 - 
 
191 
        imax=mod(mod(maxind-1,nx*ny),nx)+1
 - 
 
192 
        jmax=int(mod(maxind-1,nx*ny)/nx)+1
 - 
 
193 
        kmin=int((minind-1)/(nx*ny))+1
 - 
 
194 
        kmax=int((maxind-1)/(nx*ny))+1
 - 
 
195 
        hmin=imin+(jmin-1)*nx
 193 
      enddo
196 
        hmax=imax+(jmax-1)*ny
 - 
 
197 
        lonmin=varmin(1)+(imin-1)*dx
 194 
      lonmin = xmin + real(minindx-1) * dx
198 
        latmin=varmin(2)+(jmin-1)*dy
 195 
      latmin = ymin + real(minindy-1) * dy
199 
        pmin=ak(kmin)+bk(kmin)*sp(hmin)
 - 
 
200 
        lonmax=varmin(1)+(imax-1)*dx
 196 
      lonmax = xmin + real(maxindx-1) * dx
201 
        latmax=varmin(2)+(jmax-1)*dy
 197 
      latmax = ymin + real(maxindy-1) * dy
202 
        pmax=ak(kmax)+bk(kmax)*sp(hmax)
 - 
 
203 
      endif
 - 
 
204
198
- 
 
 199 
c     Rotate position to true longitude/latitude
205 
      if ((pollon.ne.0.).or.(pollat.ne.90.)) then
 200 
      if ((pollon.ne.0.).or.(pollat.ne.90.)) then
206 
        xphys=lmstolm(latmin,lonmin,pollat,pollon)
 201 
        xphys=lmstolm(latmin,lonmin,pollat,pollon)
207 
        yphys=phstoph(latmin,lonmin,pollat,pollon)
 202 
        yphys=phstoph(latmin,lonmin,pollat,pollon)
208 
        lonmin=xphys
 203 
        lonmin=xphys
209 
        latmin=yphys
 204 
        latmin=yphys
Line 211... Line 206...
211 
        yphys=phstoph(latmax,lonmax,pollat,pollon)
 206 
        yphys=phstoph(latmax,lonmax,pollat,pollon)
212 
        lonmax=xphys
 207 
        lonmax=xphys
213 
        latmax=yphys
 208 
        latmax=yphys
214 
      endif
 209 
      endif
215
210
216 
      enddo
 211 
c     Write output
217
- 
 
218 
      write(*,101)min,max,lonmin,latmin,nint(pmin),tmin,
 212 
      write(*,101) min,max,lonmin,latmin,nint(level),
219 
     >                    lonmax,latmax,nint(pmax),tmax
 213 
     >                     lonmax,latmax,nint(level)
220 
  101 format(2f10.3,2f8.2,i6,f6.1,2f8.2,i6,f6.1)
 214 
  101 format(2f10.3,2f8.2,i6,2f8.2,i6)
221
- 
 
222 
      call clscdf(cdfid,ierr)
 - 
 
223 
      call clscdf(cstid,ierr)
 - 
 
224
215
225 
      goto 999
 216 
c     Close netCDF file
- 
 
 217 
      call input_close(cdfid)
226
218
227 
  920 stop '*** error: variable not found on file ***'
 - 
 
228 
  999 continue
 - 
 
229 
      end
 219 
      end
230
220
231 
      subroutine vipo(var3d,varl,lev,var,nx,ny,nz,mdv,ipom)
 221 
c     ----------------------------------------------------------------------
232 
C     =====================================================
 - 
 
233
- 
 
234 
C     Interpolates the 3d variable var3d on the surface defined
 222 
C     Interpolates the 3d variable var3d on the pressure surface defined
235 
C     by lev of the variable varl. The interpolated field is
 223 
C     by lev of the variable varl. The interpolated field is
236 
C     returned as var.
 224 
C     returned as var.
237 
C     ipom determines the way of vertical interpolation:
 225 
C     ipom determines the way of vertical interpolation:
238 
C       ipom=0 is for linear interpolation
 226 
C       ipom=0 is for linear interpolation
239 
C       ipom=1 is for logarithmic interpolation
 227 
C       ipom=1 is for logarithmic interpolation
- 
 
 228 
c     ----------------------------------------------------------------------
- 
 
 229 
 
- 
 
 230 
      subroutine vipo(var3d,varl,lev,var,nx,ny,nz,mdv,ipom)
240
231 
 
241 
      integer   nx,ny,nz,ipom
 232 
      integer   nx,ny,nz,ipom
242 
      real      lev,mdv
 233 
      real      lev,mdv
243 
      real      var3d(nx,ny,nz),varl(nx,ny,nz),var(nx,ny)
 234 
      real      var3d(nx,ny,nz),varl(nx,ny,nz),var(nx,ny)
244
235
245 
      integer   i,j,k
 236 
      integer   i,j,k
Line 263... Line 254...
263 
      enddo
 254 
      enddo
264 
      enddo
 255 
      enddo
265
256
266 
      return
 257 
      return
267 
      end
 258 
      end
- 
 
 259
- 
 
 260 
c     ----------------------------------------------------------------------
- 
 
 261 
C     Interpolates the 3d variable var3d on the isentropic surface defined
- 
 
 262 
C     by lev of the variable varl. The interpolated field is
- 
 
 263 
C     returned as var.
- 
 
 264 
C     ipom determines the way of vertical interpolation:
- 
 
 265 
C       ipom=0 is for linear interpolation
- 
 
 266 
C       ipom=1 is for logarithmic interpolation
- 
 
 267 
c     ----------------------------------------------------------------------
268 
 
 268 
 
269 
      subroutine thipo(var3d,th3d,lev,var,nx,ny,nz,mdv,mode)
 269 
      subroutine thipo(var3d,th3d,lev,var,nx,ny,nz,mdv,mode)
270 
C     ======================================================
 - 
 
271
- 
 
272 
C     Interpolates the 3d variable var3d on the isentropic surface
 - 
 
273 
C     defined by lev. The interpolated field is returned as var.
 - 
 
274 
C     th3d denotes the 3d theta array.
 - 
 
275 
C     mode determines the way of vertical interpolation:
 - 
 
276 
C       mode=0 is for linear interpolation
 - 
 
277 
C       mode=1 is for logarithmic interpolation
 - 
 
278
270
279 
      integer   nx,ny,nz,mode
 271 
      integer   nx,ny,nz,mode
280 
      real      lev,mdv
 272 
      real      lev,mdv
281 
      real      var3d(nx,ny,nz),th3d(nx,ny,nz),var(nx,ny)
 273 
      real      var3d(nx,ny,nz),th3d(nx,ny,nz),var(nx,ny)
282
274
Line 301... Line 293...
301 
      enddo
 293 
      enddo
302 
      enddo
 294 
      enddo
303
295
304 
      return
 296 
      return
305 
      end
 297 
      end
- 
 
 298
- 
 
 299 
c     ----------------------------------------------------------------------
- 
 
 300 
c     Check whether character appears within string
- 
 
 301 
c     ----------------------------------------------------------------------
306 
 
 302 
 
307 
      subroutine checkchar(string,char,flag)
 303 
      subroutine checkchar(string,char,flag)
308 
C     ======================================
 - 
 
309
304
310 
      character*(*)     string
 305 
      character*(*)     string
311 
      character*(1)     char
 306 
      character*(1)     char
312 
      integer   n,flag
 307 
      integer           n,flag
313
308
314 
      flag=0
 309 
      flag=0
315 
      do n=1,len(string)
 310 
      do n=1,len(string)
316 
        if (string(n:n).eq.char) then
 311 
        if (string(n:n).eq.char) then
317 
          flag=n
 312 
          flag=n
318 
          return
 313 
          return
319 
        endif
 314 
        endif
320 
      enddo
 315 
      enddo
321 
      end
 316 
      end
- 
 
 317
- 
 
 318 
c     ----------------------------------------------------------------------
- 
 
 319 
c     3D interpolation
- 
 
 320 
c     ----------------------------------------------------------------------
322 
 
 321
323 
      real function int3d(ar,n1,n2,n3,rid,rjd,rkd)
 322 
      real function int3d(ar,n1,n2,n3,rid,rjd,rkd)
324 
c-----------------------------------------------------------------------
 - 
 
- 
 
 323 
 
325 
c     Purpose:
 324 
c     Purpose:
326 
c        This subroutine interpolates a 3d-array to an arbitrary
 325 
c        This subroutine interpolates a 3d-array to an arbitrary
327 
c        location within the grid.
 326 
c        location within the grid.
328 
c     Arguments:
 327 
c     Arguments:
329 
c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
 328 
c        ar        real  input   surface pressure, define as ar(n1,n2,n3)
330 
c        n1,n2,n3  int   input   dimensions of ar
 329 
c        n1,n2,n3  int   input   dimensions of ar
331 
c        ri,rj,rk  real  input   grid location to be interpolated to
 330 
c        ri,rj,rk  real  input   grid location to be interpolated to
332 
c     History:
 331 
c     History:
333 
c-----------------------------------------------------------------------
 - 
 
334 
 
 332 
 
335 
c     argument declarations
 333 
c     argument declarations
336 
      integer   n1,n2,n3
 334 
      integer   n1,n2,n3
337 
      real      ar(n1,n2,n3), rid,rjd,rkd
 335 
      real      ar(n1,n2,n3), rid,rjd,rkd
338 
 
 336 
 
Line 417... Line 415...
417 
     &           + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
 415 
     &           + ar(ip1,jp1,kp1) * frac0i * frac0j * frac0k
418 
c          print *,'int3d 11: ',rid,rjd,rkd,int3d
 416 
c          print *,'int3d 11: ',rid,rjd,rkd,int3d
419 
        endif
 417 
        endif
420 
      endif
 418 
      endif
421 
      end
 419 
      end
- 
 
 420
- 
 
 421 
c     ----------------------------------------------------------------------
- 
 
 422 
c     3D interpolation including missing data check
- 
 
 423 
c     ----------------------------------------------------------------------
- 
 
 424
422 
      real function int3dm(ar,n1,n2,n3,rid,rjd,rkd,misdat)
 425 
      real function int3dm(ar,n1,n2,n3,rid,rjd,rkd,misdat)
423 
c-----------------------------------------------------------------------
 - 
 
- 
 
 426 
 
424 
c     Purpose:
 427 
c     Purpose:
425 
c        This subroutine interpolates a 3d-array to an arbitrary
 428 
c        This subroutine interpolates a 3d-array to an arbitrary
426 
c        location within the grid. The interpolation includes the
429 
c        location within the grid. The interpolation includes the
427 
c        testing of the missing data flag 'misdat'.
430 
c        testing of the missing data flag 'misdat'.
428 
c     Arguments:
 431 
c     Arguments:
Line 431... Line 434...
431 
c        ri,rj,rk  real  input   grid location to be interpolated to
 434 
c        ri,rj,rk  real  input   grid location to be interpolated to
432 
c        misdat    real  input   missing data flag (on if misdat<>0)
 435 
c        misdat    real  input   missing data flag (on if misdat<>0)
433 
c     Warning:
 436 
c     Warning:
434 
c        This routine has not yet been seriously tested
 437 
c        This routine has not yet been seriously tested
435 
c     History:
 438 
c     History:
436 
c-----------------------------------------------------------------------
 - 
 
437 
 
 439 
 
438 
c     argument declarations
 440 
c     argument declarations
439 
      integer   n1,n2,n3
 441 
      integer   n1,n2,n3
440 
      real      ar(n1,n2,n3), rid,rjd,rkd, misdat
 442 
      real      ar(n1,n2,n3), rid,rjd,rkd, misdat
441 
 
 443 
 
Line 554... Line 556...
554 
           endif
 556 
           endif
555 
        endif
 557 
        endif
556 
      endif
 558 
      endif
557 
      end
 559 
      end
558 
 
 560 
 
- 
 
 561 
c     ----------------------------------------------------------------------
- 
 
 562 
c     Calculate potential temperature
- 
 
 563 
c     ----------------------------------------------------------------------
559 
 
 564 
 
560 
      subroutine pottemp(pt,t,sp,ie,je,ke,ak,bk)
 565 
      subroutine pottemp(pt,t,sp,ie,je,ke,ak,bk)
561 
c     ==========================================
 - 
 
562
566
563 
c     argument declaration
 567 
c     argument declaration
564 
      integer   ie,je,ke
 568 
      integer   ie,je,ke
565 
      real      pt(ie,je,ke),t(ie,je,ke),sp(ie,je),
 569 
      real      pt(ie,je,ke),t(ie,je,ke),sp(ie,je),
566 
     >     	ak(ke),bk(ke)
 570 
     >     	ak(ke),bk(ke)
Line 588... Line 592...
588 
          endif
 592 
          endif
589 
        enddo
 593 
        enddo
590 
      enddo
 594 
      enddo
591 
      enddo
 595 
      enddo
592 
      end
 596 
      end
- 
 
 597
- 
 
 598 
c     ----------------------------------------------------------------------
- 
 
 599 
c     Calculate 3D pressure
- 
 
 600 
c     ----------------------------------------------------------------------
593 
 
 601 
 
594 
      subroutine pres(pr,sp,ie,je,ke,ak,bk)
 602 
      subroutine pres(pr,sp,ie,je,ke,ak,bk)
595 
c     =====================================
 - 
 
- 
 
 603 
 
596 
c     argument declaration
 604 
c     argument declaration
597 
      integer  ie,je,ke
 605 
      integer  ie,je,ke
598 
      real,intent(OUT) :: pr(ie,je,ke)
 606 
      real,intent(OUT) :: pr(ie,je,ke)
599 
      real,intent(IN)  :: sp(ie,je)
 607 
      real,intent(IN)  :: sp(ie,je)
600 
      real,intent(IN)  :: ak(ke),bk(ke)
 608 
      real,intent(IN)  :: ak(ke),bk(ke)
Line 609... Line 617...
609 
        pr(i,j,k)=ak(k)+bk(k)*sp(i,j)
 617 
        pr(i,j,k)=ak(k)+bk(k)*sp(i,j)
610 
      enddo
 618 
      enddo
611 
      enddo
 619 
      enddo
612 
      enddo
 620 
      enddo
613 
      end
 621 
      end
- 
 
 622
- 
 
 623 
c     ----------------------------------------------------------------------
- 
 
 624 
c     Coordinate rotations from COSMO
- 
 
 625 
c     ----------------------------------------------------------------------
- 
 
 626
614 
      REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
 627 
      REAL FUNCTION PHTOPHS (PHI, LAM, POLPHI, POLLAM)
615 
C
 628 
C
616 
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
 629 
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
617 
C
 630 
C
618 
C**** PHTOPHS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
 631 
C**** PHTOPHS  -   FC:UMRECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE PHI
Line 660... Line 673...
660
673
661 
      PHTOPHS = ZRPI18*ASIN(ZARG)
 674 
      PHTOPHS = ZRPI18*ASIN(ZARG)
662
675
663 
      RETURN
 676 
      RETURN
664 
      END
 677 
      END
- 
 
 678
- 
 
 679
665 
      REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
 680 
      REAL FUNCTION PHSTOPH (PHIS, LAMS, POLPHI, POLLAM)
666 
C
 681 
C
667 
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
 682 
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
668 
C
 683 
C
669 
C**** PHSTOPH  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
 684 
C**** PHSTOPH  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN BREITE FUER
Line 772... Line 787...
772 
        LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
 787 
        LMTOLMS = ZRPI18*ATAN2(ZARG1,ZARG2)
773 
      ENDIF
 788 
      ENDIF
774
789
775 
      RETURN
 790 
      RETURN
776 
      END
 791 
      END
- 
 
 792
- 
 
 793
777 
      REAL FUNCTION LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
 794 
      REAL FUNCTION LMSTOLM (PHIS, LAMS, POLPHI, POLLAM)
778 
C
 795 
C
779 
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
 796 
C%Z% Modul %M%, V%I% vom %G%, extrahiert am %H%
780 
C
 797 
C
781 
C**** LMSTOLM  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER
 798 
C**** LMSTOLM  -   FC:BERECHNUNG DER WAHREN GEOGRAPHISCHEN LAENGE FUER