Subversion Repositories lagranto.ecmwf

Rev

Rev 13 | Rev 23 | Go to most recent revision | Show entire file | Regard 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 
 
 8 
 
11 
      implicit none
 9 
      use netcdf
12 
 
 10 
 
13 
      integer	nzmax,ntmax
 11 
      implicit none
14 
      parameter(nzmax=100,ntmax=200)
 - 
 
15 
 
 12 
 
- 
 
 13 
c     ----------------------------------------------------------------------
- 
 
 14 
c     Declaration of parameters and variables
- 
 
 15 
c     ----------------------------------------------------------------------
- 
 
 16
- 
 
 17 
c     Interpolation method
- 
 
 18 
      integer   ipom
- 
 
 19 
      parameter (ipom=0)
- 
 
 20
- 
 
 21 
c     Flag for timecheck
- 
 
 22 
      character*80 timecheck
- 
 
 23 
      parameter    (timecheck = 'no' )
- 
 
 24
- 
 
 25 
c     netCDF fields
- 
 
 26 
      real,allocatable, dimension (:,:)   :: sp,varf
16 
      real,allocatable, dimension (:) :: sp,varf,field,varl,tt
 27 
      real,allocatable, dimension (:,:,:) :: field,varl,tt,p3
- 
 
 28 
      real,allocatable, dimension (:)     :: ak,bk
17 
      integer   stat
 29 
      integer                                stat
- 
 
 30 
      character*80                           cdfnam,varnam
- 
 
 31 
      character*80                           vnam(200)
- 
 
 32 
      integer                                nvars
- 
 
 33 
      integer	                             cdfid,ierr
- 
 
 34 
      real	                                 mdv
- 
 
 35 
      real                                   stagz
- 
 
 36
- 
 
 37 
c     Grid description
- 
 
 38 
      real      xmin,xmax       ! Zonal grid extension
- 
 
 39 
      real      ymin,ymax       ! Meridional grid extension
- 
 
 40 
      integer   nx,ny,nz        ! Grid dimensions
- 
 
 41 
      real      dx,dy           ! Horizontal grid resolution
- 
 
 42 
      real	    pollon,pollat   ! Pole position
- 
 
 43 
 
- 
 
 44 
c     Output variables
- 
 
 45 
      real      min,max                      ! Minimum & maximum value
- 
 
 46 
      real      lonmin,lonmax,latmin,latmax  ! Position of min & max
18 
 
 47 
 
- 
 
 48 
c     Auxiliary variables
- 
 
 49 
      integer        iargc
- 
 
 50 
      character*(80) arg
19 
      real	time(ntmax),level
 51 
      real           rd
- 
 
 52 
      integer        i,j
20 
      character*80 cdfnam,cstnam,varnam,tvar,clev
 53 
      integer        minindx,minindy,maxindx,maxindy
21 
      character*20 vnam(100)
 54 
      character*80   tvar
22 
      character*1  mode
 55 
      character*1    mode
23 
      integer	cdfid,cstid,ierr,ndim,vardim(4)
 - 
 
24 
      real	dx,dy,mdv,varmin(4),varmax(4),stag(4)
 - 
 
25 
      real      aklev(nzmax),bklev(nzmax),aklay(nzmax),bklay(nzmax),
 - 
 
26 
     >          ak(nzmax),bk(nzmax)
 - 
 
27 
      logical	prelev
 56 
      character*80   clev
28 
      integer	nx,ny,nz,ntimes,ipom,nvars
 57 
      real	         xphys,yphys
29 
      real	pollon,pollat,xphys,yphys
 58 
      real	         level
30 
      integer	i,n
 59 
      integer        flag
31 
      integer	imin,jmin,imax,jmax,minind,maxind,kmin,kmax,hmin,hmax
 - 
 
32 
      real	min,max,lonmin,latmin,lonmax,latmax,pmin,pmax,tmin,tmax
 - 
 
33 
 
 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)
 - 
 
74 
      call getcfn(cdfid,cstnam,ierr)
 - 
 
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
- 
 
80 
C     Get the data array and the levels
 101 
c     Get list of variables on file
81
- 
 
82 
      nx=vardim(1)
 - 
 
83 
      ny=vardim(2)
 - 
 
84 
      nz=vardim(3)
 - 
 
85 
      call gettimes(cdfid,time,ntimes,ierr)
 - 
 
86 
 
 - 
 
87 
      call getgrid(cstid,dx,dy,ierr)
 102 
      call input_getvars (cdfid,vnam,nvars)
88 
      call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr)
 - 
 
89 
      call getpole(cstid,pollon,pollat,ierr)
 - 
 
90 
 
 103 
 
91 
C     Get memory for dynamic arrays
 104 
C     Get infos about data domain - in particular: dimensions
- 
 
 105 
      nx       = 1
- 
 
 106 
      ny       = 1
- 
 
 107 
      nz       = 1
- 
 
 108 
      call input_grid (-cdfid,varnam,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
- 
 
 109 
     >               0.,pollon,pollat,rd,rd,nz,rd,rd,rd,stagz,timecheck)
92
110 
 
- 
 
 111 
C     Get memory for dynamic arrays
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 ***'
- 
 
 126
- 
 
 127 
c     ----------------------------------------------------------------------
- 
 
 128 
c     Load data
- 
 
 129 
c     ----------------------------------------------------------------------
- 
 
 130 
 
- 
 
 131 
c     Get grid info for variable - non-constant part  (P, PS, AK, BK)
- 
 
 132 
      call input_grid
- 
 
 133 
     >       (cdfid,varnam,xmin,xmax,ymin,ymax,dx,dy,nx,ny,
- 
 
 134 
     >        0.,pollon,pollat,varl,sp,nz,ak,bk,stagz,timecheck)
- 
 
 135
- 
 
 136 
c     Load Variable
- 
 
 137 
      call input_wind (cdfid,varnam,field,0.,stagz,mdv,
- 
 
 138 
     >                 xmin,xmax,ymin,ymax,dx,dy,nx,ny,nz,timecheck)
103
139
104 
      do n=1,ntimes
 - 
 
105
- 
 
106 
      call getdat(cdfid,trim(varnam),time(n),0,field,ierr)
 - 
 
107
- 
 
108 
C     Define the appropriate ak and bk-arrays
 140 
c     In case of isentropic level, read temperature
109
- 
 
110 
      if (stag(3).eq.0.) then
 141 
      if (mode.eq.'T') then
111 
        do i=1,nz
 - 
 
112 
          ak(i)=aklev(i)
 - 
 
113 
          bk(i)=bklev(i)
 - 
 
114 
        enddo
 - 
 
115 
      else
 - 
 
116 
        do i=1,nz
 - 
 
117 
          ak(i)=aklay(i)
 - 
 
118 
          bk(i)=bklay(i)
 - 
 
119 
        enddo
 - 
 
120 
      endif
 - 
 
121 
 
 - 
 
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)
 - 
 
127 
 
 - 
 
128 
C     Determine name of "temperature variable"
 - 
 
129 
 
 - 
 
130 
      call getvars(cdfid,nvars,vnam,ierr)
 - 
 
131 
      tvar="T"
 142 
            tvar="T"
132 
      do i=1,nvars
 143 
            do i=1,nvars
133 
        if (vnam(i).eq."THETA") tvar="THETA"
 144 
              if (vnam(i).eq."THETA") tvar="THETA"
134 
        if (vnam(i).eq."TH")    tvar="TH"
 145 
              if (vnam(i).eq."TH")    tvar="TH"
135 
      enddo
 146 
            enddo
- 
 
 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)
- 
 
 151 
            else
- 
 
 152 
                 call input_wind (cdfid,tvar,tt,0.,stagz,mdv,xmin,xmax,
- 
 
 153 
     >                            ymin,ymax,dx,dy,nx,ny,nz,timecheck)
- 
 
 154 
            endif
- 
 
 155 
      endif
136 
 
 156 
 
- 
 
 157 
c     ----------------------------------------------------------------------
137 
C     If required do interpolation on pressure or theta level
 158 
c     Do interpolation and get minimum, maximum
- 
 
 159 
c     ----------------------------------------------------------------------
138
160 
 
- 
 
 161 
C     If required do interpolation on pressure or theta level
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
143 
        if (tvar.eq.'T') then
 - 
 
144 
          call getdat(cdfid,'T',time(n),0,tt,ierr)
 - 
 
145 
          call pottemp(varl,tt,sp,nx,ny,nz,ak,bk)
 - 
 
146 
          call thipo(field,varl,level,varf,nx,ny,nz,mdv,ipom)
 165 
        call thipo(field,varl,level,varf,nx,ny,nz,mdv,ipom)
- 
 
 166 
      elseif (mode.eq.'L') then
147 
        else
 167 
        do i=1,nx
- 
 
 168 
           do j=1,ny
148 
          call getdat(cdfid,tvar,time(n),0,varl,ierr)
 169 
               varf(i,j) = field(i,j,nint(level))
- 
 
 170 
           enddo
- 
 
 171 
        enddo
- 
 
 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
 - 
 
188 
        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
 - 
 
196 
        hmax=imax+(jmax-1)*ny
 - 
 
197 
        lonmin=varmin(1)+(imin-1)*dx
 - 
 
198 
        latmin=varmin(2)+(jmin-1)*dy
 - 
 
199 
        pmin=ak(kmin)+bk(kmin)*sp(hmin)
 - 
 
200 
        lonmax=varmin(1)+(imax-1)*dx
 - 
 
201 
        latmax=varmin(2)+(jmax-1)*dy
 - 
 
202 
        pmax=ak(kmax)+bk(kmax)*sp(hmax)
 - 
 
203 
      endif
 191 
          endif
- 
 
 192 
        enddo
- 
 
 193 
      enddo
- 
 
 194 
      lonmin = xmin + real(minindx-1) * dx
- 
 
 195 
      latmin = ymin + real(minindy-1) * dy
- 
 
 196 
      lonmax = xmin + real(maxindx-1) * dx
- 
 
 197 
      latmax = ymin + real(maxindy-1) * dy
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
215
222 
      call clscdf(cdfid,ierr)
 216 
c     Close netCDF file
223 
      call clscdf(cstid,ierr)
 217 
      call input_close(cdfid)
224
218
225 
      goto 999
 - 
 
226
- 
 
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)
 - 
 
232 
C     =====================================================
 221 
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
Line 264... Line 255...
264 
      enddo
 255 
      enddo
265
256
266 
      return
 257 
      return
267 
      end
 258 
      end
268 
 
 259
- 
 
 260 
c     ----------------------------------------------------------------------
269 
      subroutine thipo(var3d,th3d,lev,var,nx,ny,nz,mdv,mode)
 261 
C     Interpolates the 3d variable var3d on the isentropic surface defined
270 
C     ======================================================
 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     ----------------------------------------------------------------------
271
268 
 
272 
C     Interpolates the 3d variable var3d on the isentropic surface
 269 
      subroutine thipo(var3d,th3d,lev,var,nx,ny,nz,mdv,mode)
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 302... Line 294...
302 
      enddo
 294 
      enddo
303
295
304 
      return
 296 
      return
305 
      end
 297 
      end
306 
 
 298
- 
 
 299 
c     ----------------------------------------------------------------------
- 
 
 300 
c     Check whether character appears within string
- 
 
 301 
c     ----------------------------------------------------------------------
- 
 
 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
Line 318... Line 313...
318 
          return
 313 
          return
319 
        endif
 314 
        endif
320 
      enddo
 315 
      enddo
321 
      end
 316 
      end
322 
 
 317
- 
 
 318 
c     ----------------------------------------------------------------------
- 
 
 319 
c     3D interpolation
- 
 
 320 
c     ----------------------------------------------------------------------
- 
 
 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 589... Line 593...
589 
        enddo
 593 
        enddo
590 
      enddo
 594 
      enddo
591 
      enddo
 595 
      enddo
592 
      end
 596 
      end
593 
 
 597
- 
 
 598 
c     ----------------------------------------------------------------------
- 
 
 599 
c     Calculate 3D pressure
- 
 
 600 
c     ----------------------------------------------------------------------
- 
 
 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