getfields.f90 9.73 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
151
152
153
        nstop = 1
        goto 40
      endif
Matthias Langer's avatar
 
Matthias Langer committed
154
    end do
155
40  indmin=indj
Matthias Langer's avatar
 
Matthias Langer committed
156

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

Matthias Langer's avatar
 
Matthias Langer committed
161
162
  else

163
164
165
166
167
168
169
170
! 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
171
172
173
174
175
        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
176
        call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn)
177
178
179
180
181
182
183
        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
184
185
186
        call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
        memtime(1)=wftime(indj)
        memind(2)=2
187
188
189
190
191
        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
192
        call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
193
194
195
196
197
198
199
        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
200
201
202
203
204
205
206
        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
207

208
209
210
211
   if (WETBKDEP) then
        call writeprecip(itime,memind(1))
   endif

Matthias Langer's avatar
 
Matthias Langer committed
212
213
214
215
216
217
218
  endif

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

  if (lwindinterv.gt.idiffmax) nstop=3

end subroutine getfields