getfields.f90 10.1 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
        call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
        memtime(2)=wftime(indj+1)
149
        if (WATERCYCLE) &
150
           call euler_rain(itime) ! calculate the integrated E-P
Sabine's avatar
Sabine committed
151
           call euler_rain_nest(itime) ! calculate the integrated E-P
152
153
154
        nstop = 1
        goto 40
      endif
Matthias Langer's avatar
 
Matthias Langer committed
155
    end do
156
40  indmin=indj
Matthias Langer's avatar
 
Matthias Langer committed
157

158
   if ((WETBKDEP).or.(WATERCYCLE)) then
159
160
161
        call writeprecip(itime,memind(1))
   endif

Matthias Langer's avatar
 
Matthias Langer committed
162
163
  else

164
165
166
167
168
169
170
171
! 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
172
173
174
175
176
        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
177
        call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn)
178
179
180
181
182
183
184
        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
185
186
187
        call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
        memtime(1)=wftime(indj)
        memind(2)=2
Sabine's avatar
Sabine committed
188
        if (WATERCYCLE) then
189
           call euler_rain(itime) ! calculate the integrated E-P
Sabine's avatar
Sabine committed
190
191
           call euler_rain_nest(itime) ! calculate the integrated E-P
        endif
192
193
194
195
196
        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
197
        call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
198
199
200
201
202
203
204
        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
205
206
        call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
        memtime(2)=wftime(indj+1)
Sabine's avatar
Sabine committed
207
        if (WATERCYCLE) then
208
           call euler_rain(itime) ! calculate the integrated E-P
Sabine's avatar
Sabine committed
209
210
           call euler_rain_nest(itime) ! calculate the integrated E-P
        endif
211
212
213
214
215
        nstop = 1
        goto 60
      endif
    end do
60  indmin=indj
Matthias Langer's avatar
 
Matthias Langer committed
216

217
218
219
220
   if (WETBKDEP) then
        call writeprecip(itime,memind(1))
   endif

Matthias Langer's avatar
 
Matthias Langer committed
221
222
223
224
225
226
227
  endif

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

  if (lwindinterv.gt.idiffmax) nstop=3

end subroutine getfields