getfields.f90 9.62 KB
Newer Older
Matthias Langer's avatar
   
Matthias Langer committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
!**********************************************************************
! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
!                                                                     *
! This file is part of FLEXPART.                                      *
!                                                                     *
! FLEXPART is free software: you can redistribute it and/or modify    *
! it under the terms of the GNU General Public License as published by*
! the Free Software Foundation, either version 3 of the License, or   *
! (at your option) any later version.                                 *
!                                                                     *
! FLEXPART is distributed in the hope that it will be useful,         *
! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
! GNU General Public License for more details.                        *
!                                                                     *
! You should have received a copy of the GNU General Public License   *
! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
!**********************************************************************

22
subroutine getfields(itime,nstop,metdata_format)
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
!                       i     o
!*****************************************************************************
!                                                                            *
!  This subroutine manages the 3 data fields to be kept in memory.           *
!  During the first time step of petterssen it has to be fulfilled that the  *
!  first data field must have |wftime|<itime, i.e. the absolute value of     *
!  wftime must be smaller than the absolute value of the current time in [s].*
!  The other 2 fields are the next in time after the first one.              *
!  Pointers (memind) are used, because otherwise one would have to resort the*
!  wind fields, which costs a lot of computing time. Here only the pointers  *
!  are resorted.                                                             *
!                                                                            *
!     Author: A. Stohl                                                       *
!                                                                            *
!     29 April 1994                                                          *
!                                                                            *
!*****************************************************************************
40
41
42
43
44
45
46
!  Changes, Bernd C. Krueger, Feb. 2001:                                     *
!   Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block.        *
!   Function of nstop extended.                                              *
!                                                                            *
!   Unified ECMWF and GFS builds                                             *
!   Marian Harustak, 12.5.2017                                               *
!     - Added passing of metdata_format as it was needed by called routines  *
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
!*****************************************************************************
!                                                                            *
! Variables:                                                                 *
! lwindinterval [s]    time difference between the two wind fields read in   *
! indj                 indicates the number of the wind field to be read in  *
! indmin               remembers the number of wind fields already treated   *
! memind(2)            pointer, on which place the wind fields are stored    *
! memtime(2) [s]       times of the wind fields, which are kept in memory    *
! itime [s]            current time since start date of trajectory calcu-    *
!                      lation                                                *
! nstop                > 0, if trajectory has to be terminated               *
! nx,ny,nuvz,nwz       field dimensions in x,y and z direction               *
! uu(0:nxmax,0:nymax,nuvzmax,2)   wind components in x-direction [m/s]       *
! vv(0:nxmax,0:nymax,nuvzmax,2)   wind components in y-direction [m/s]       *
! ww(0:nxmax,0:nymax,nwzmax,2)    wind components in z-direction [deltaeta/s]*
! tt(0:nxmax,0:nymax,nuvzmax,2)   temperature [K]                            *
! ps(0:nxmax,0:nymax,2)           surface pressure [Pa]                      *
64
! metdata_format     format of metdata (ecmwf/gfs)                           *
65
66
67
68
69
70
!                                                                            *
! Constants:                                                                 *
! idiffmax             maximum allowable time difference between 2 wind      *
!                      fields                                                *
!                                                                            *
!*****************************************************************************
Matthias Langer's avatar
   
Matthias Langer committed
71
72
73

  use par_mod
  use com_mod
74
  use class_gribfile
Matthias Langer's avatar
   
Matthias Langer committed
75
76
77
78

  implicit none

  integer :: indj,itime,nstop,memaux
79
  integer :: metdata_format
Matthias Langer's avatar
   
Matthias Langer committed
80
81
82
83
84
85
86
87
88
89
90
91
92

  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
  real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)

  integer :: indmin = 1


93
94
! Check, if wind fields are available for the current time step
!**************************************************************
Matthias Langer's avatar
   
Matthias Langer committed
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109

  nstop=0

  if ((ldirect*wftime(1).gt.ldirect*itime).or. &
       (ldirect*wftime(numbwf).lt.ldirect*itime)) then
    write(*,*) 'FLEXPART WARNING: NO WIND FIELDS ARE AVAILABLE.'
    write(*,*) 'A TRAJECTORY HAS TO BE TERMINATED.'
    nstop=4
    return
  endif


  if ((ldirect*memtime(1).le.ldirect*itime).and. &
       (ldirect*memtime(2).gt.ldirect*itime)) then

110
111
! The right wind fields are already in memory -> don't do anything
!*****************************************************************
Matthias Langer's avatar
   
Matthias Langer committed
112
113
114
115
116
117
118

    continue

  else if ((ldirect*memtime(2).le.ldirect*itime).and. &
       (memtime(2).ne.999999999)) then


119
120
121
! Current time is after 2nd wind field
! -> Resort wind field pointers, so that current time is between 1st and 2nd
!***************************************************************************
Matthias Langer's avatar
   
Matthias Langer committed
122
123
124
125
126
127
128

    memaux=memind(1)
    memind(1)=memind(2)
    memind(2)=memaux
    memtime(1)=memtime(2)


129
130
! Read a new wind field and store it on place memind(2)
!******************************************************
Matthias Langer's avatar
   
Matthias Langer committed
131
132

    do indj=indmin,numbwf-1
133
      if (ldirect*wftime(indj+1).gt.ldirect*itime) then
134
135
136
137
138
        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
          call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
        else
          call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh)
        end if
139
        call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
140
141
142
143
144
145
146
        call calcpar(memind(2),uuh,vvh,pvh,metdata_format)
        call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format)
        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
          call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh)
        else
          call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh)
        end if
147
148
149
150
151
        call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
        memtime(2)=wftime(indj+1)
        nstop = 1
        goto 40
      endif
Matthias Langer's avatar
   
Matthias Langer committed
152
    end do
153
40  indmin=indj
Matthias Langer's avatar
   
Matthias Langer committed
154

155
156
157
158
   if (WETBKDEP) then
        call writeprecip(itime,memind(1))
   endif

Matthias Langer's avatar
   
Matthias Langer committed
159
160
  else

161
162
163
164
165
166
167
168
! No wind fields, which can be used, are currently in memory
! -> read both wind fields
!***********************************************************

    do indj=indmin,numbwf-1
      if ((ldirect*wftime(indj).le.ldirect*itime).and. &
           (ldirect*wftime(indj+1).gt.ldirect*itime)) then
        memind(1)=1
169
170
171
172
173
        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
          call readwind_ecmwf(indj,memind(1),uuh,vvh,wwh)
        else
          call readwind_gfs(indj,memind(1),uuh,vvh,wwh)
        end if
174
        call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn)
175
176
177
178
179
180
181
        call calcpar(memind(1),uuh,vvh,pvh,metdata_format)
        call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format)
        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
          call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh)
        else
          call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh)
        end if
182
183
184
        call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
        memtime(1)=wftime(indj)
        memind(2)=2
185
186
187
188
189
        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
          call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
        else
          call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh)
        end if
190
        call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
191
192
193
194
195
196
197
        call calcpar(memind(2),uuh,vvh,pvh,metdata_format)
        call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format)
        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
          call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh)
        else
          call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh)
        end if
198
199
200
201
202
203
204
        call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
        memtime(2)=wftime(indj+1)
        nstop = 1
        goto 60
      endif
    end do
60  indmin=indj
Matthias Langer's avatar
   
Matthias Langer committed
205

206
207
208
209
   if (WETBKDEP) then
        call writeprecip(itime,memind(1))
   endif

Matthias Langer's avatar
   
Matthias Langer committed
210
211
212
213
214
215
216
  endif

  lwindinterv=abs(memtime(2)-memtime(1))

  if (lwindinterv.gt.idiffmax) nstop=3

end subroutine getfields