source: flexpart.git/src/getfields_mpi.f90 @ 02095e3

devrelease-10univie
Last change on this file since 02095e3 was 6ecb30a, checked in by Espen Sollum ATMOS <eso@…>, 2 years ago

Merged changes from CTBTO project

  • Property mode set to 100644
File size: 12.1 KB
Line 
1!**********************************************************************
2! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
3! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
4! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
5!                                                                     *
6! This file is part of FLEXPART.                                      *
7!                                                                     *
8! FLEXPART is free software: you can redistribute it and/or modify    *
9! it under the terms of the GNU General Public License as published by*
10! the Free Software Foundation, either version 3 of the License, or   *
11! (at your option) any later version.                                 *
12!                                                                     *
13! FLEXPART is distributed in the hope that it will be useful,         *
14! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
15! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
16! GNU General Public License for more details.                        *
17!                                                                     *
18! You should have received a copy of the GNU General Public License   *
19! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
20!**********************************************************************
21
22subroutine getfields(itime,nstop,memstat,metdata_format)
23!                       i     o       o
24!*****************************************************************************
25!                                                                            *
26!  This subroutine manages the 3 data fields to be kept in memory.           *
27!  During the first time step of petterssen it has to be fulfilled that the  *
28!  first data field must have |wftime|<itime, i.e. the absolute value of     *
29!  wftime must be smaller than the absolute value of the current time in [s].*
30!  The other 2 fields are the next in time after the first one.              *
31!  Pointers (memind) are used, because otherwise one would have to resort the*
32!  wind fields, which costs a lot of computing time. Here only the pointers  *
33!  are resorted.                                                             *
34!                                                                            *
35!     Author: A. Stohl                                                       *
36!                                                                            *
37!     29 April 1994                                                          *
38!                                                                            *
39!*****************************************************************************
40!  Changes, Bernd C. Krueger, Feb. 2001:
41!   Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block.
42!   Function of nstop extended.
43!
44!  eso 2014:
45!    MPI version.
46!    If running with number of processes >= mpi_mod::read_grp_min,
47!    only one process (mpi_mod::lmpreader=.true.) does the actual reading, while
48!    the rest call this routine only to update memind, memstat etc.
49!
50!    If mpi_mod::lmp_sync=.true., uses 3 fields instead of 2, to allow reading
51!    the newest in the background.
52!
53!    Return memstat, which is the sum of
54!   
55!    memstat=16:     one new field to be read
56!    memstat=32:     two new fields must be read
57!    memstat=1,2,3:  position(s) in array to read next field
58!    memstat=0:      no new fields to be read
59!         
60!   Unified ECMWF and GFS builds                                             
61!   Marian Harustak, 12.5.2017                                               
62!     - Added passing of metdata_format as it was needed by called routines 
63!         
64!*****************************************************************************
65!                                                                            *
66! Variables:                                                                 *
67! lwindinterval [s]    time difference between the two wind fields read in   *
68! indj                 indicates the number of the wind field to be read in  *
69! indmin               remembers the number of wind fields already treated   *
70! memind(2[3])         pointer, on which place the wind fields are stored    *
71! memtime(2[3]) [s]    times of the wind fields, which are kept in memory    *
72! itime [s]            current time since start date of trajectory calcu-    *
73!                      lation                                                *
74! nstop                > 0, if trajectory has to be terminated               *
75! memstat              keep track of change of status for field pointers     *
76! nx,ny,nuvz,nwz       field dimensions in x,y and z direction               *
77! uu(0:nxmax,0:nymax,nuvzmax,numwfmem)   wind components in x-direction [m/s]*
78! vv(0:nxmax,0:nymax,nuvzmax,numwfmem)   wind components in y-direction [m/s]*
79! ww(0:nxmax,0:nymax,nwzmax,numwfmem)    wind components in z-direction      *
80!                                          [deltaeta/s]                      *
81! tt(0:nxmax,0:nymax,nuvzmax,numwfmem)   temperature [K]                     *
82! ps(0:nxmax,0:nymax,numwfmem)           surface pressure [Pa]               *
83! metdata_format     format of metdata (ecmwf/gfs)                           *
84!                                                                            *
85! Constants:                                                                 *
86! idiffmax             maximum allowable time difference between 2 wind      *
87!                      fields                                                *
88!                                                                            *
89!*****************************************************************************
90
91  use par_mod
92  use com_mod
93  use mpi_mod, only: lmpreader,lmp_use_reader,lmp_sync
94  use class_gribfile
95
96  implicit none
97
98  integer :: metdata_format
99  integer :: indj,itime,nstop,memaux,mindread
100! mindread: index of where to read 3rd field
101  integer, intent(out) :: memstat
102! keep track of 3rd field index. updated when new fields are read
103  integer, save :: mind3=0
104
105  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
106  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
107  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
108  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
109  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
110  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
111  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
112  real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
113
114  integer :: indmin = 1
115
116
117! Check, if wind fields are available for the current time step
118!**************************************************************
119
120  nstop=0
121  memstat=0
122
123  if ((ldirect*wftime(1).gt.ldirect*itime).or. &
124       (ldirect*wftime(numbwf).lt.ldirect*itime)) then
125    write(*,*) 'FLEXPART WARNING: NO WIND FIELDS ARE AVAILABLE.'
126    write(*,*) 'A TRAJECTORY HAS TO BE TERMINATED.'
127    nstop=4
128    return
129  endif
130
131
132  if ((ldirect*memtime(1).le.ldirect*itime).and. &
133       (ldirect*memtime(2).gt.ldirect*itime)) then
134
135! The right wind fields are already in memory -> don't do anything
136!*****************************************************************
137    memstat=0
138    continue
139
140  else if ((ldirect*memtime(2).le.ldirect*itime).and. &
141       (memtime(2).ne.999999999)) then
142
143
144! Current time is after 2nd wind field
145! -> Resort wind field pointers, so that current time is between 1st and 2nd
146!***************************************************************************
147
148! 2 fields in memory
149!*******************
150    if (lmp_sync) then
151      memaux=memind(1)
152      memind(1)=memind(2)
153      memind(2)=memaux
154      memtime(1)=memtime(2)
155      memstat=16+memind(2)
156      mindread=memind(2)
157
158    else
159! 3 fields in memory
160! Note that the reader process never needs to keep 3 fields in memory,
161! (2 is enough) so it is possible to save some memory
162!*********************************************************************
163      if (mind3.eq.32.or.mind3.eq.19) then
164        if (lmpreader) then
165          memind(1)=2
166          memind(2)=3
167          memind(3)=3
168        else
169          memind(1)=2
170          memind(2)=3
171          memind(3)=1
172        end if
173        memstat=17
174      else if (mind3.eq.17) then
175        if (lmpreader) then
176          memind(1)=3
177          memind(2)=1
178          memind(3)=1
179        else
180          memind(1)=3
181          memind(2)=1
182          memind(3)=2
183        end if
184        memstat=18
185      else if (mind3.eq.18) then
186        if (lmpreader) then
187          memind(1)=1
188          memind(2)=2
189          memind(3)=2
190        else
191          memind(1)=1
192          memind(2)=2
193          memind(3)=3
194        end if
195        memstat=19
196      else
197        write(*,*) '#### getfields_mpi.f90> ERROR: ', &
198             & 'unknown mind3, exiting.', mind3,' ####'
199        stop
200      end if
201      mind3=memstat
202      memtime(1)=memtime(2)
203      mindread=memind(3)
204    end if
205
206
207! Read a new wind field and store it on place memind(2)
208! or memind(3) for the 3-field read-ahead version
209!******************************************************
210    do indj=indmin,numbwf-1
211      if (ldirect*wftime(indj+1).gt.ldirect*itime) then
212        if (lmpreader.or..not. lmp_use_reader) then
213          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
214            call readwind_ecmwf(indj+1,mindread,uuh,vvh,wwh)
215          else
216            call readwind_gfs(indj+1,mindread,uuh,vvh,wwh)
217          end if
218          call readwind_nests(indj+1,mindread,uuhn,vvhn,wwhn)
219          call calcpar(mindread,uuh,vvh,pvh,metdata_format)
220          call calcpar_nests(mindread,uuhn,vvhn,pvhn,metdata_format)
221          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
222            call verttransform_ecmwf(mindread,uuh,vvh,wwh,pvh)
223          else
224            call verttransform_gfs(mindread,uuh,vvh,wwh,pvh)
225          end if
226          call verttransform_nests(mindread,uuhn,vvhn,wwhn,pvhn)
227        end if
228        memtime(2)=wftime(indj+1)
229        nstop = 1
230        goto 40
231      endif
232    end do
23340  indmin=indj
234
235  else
236
237! No wind fields, which can be used, are currently in memory
238! -> read both wind fields
239!***********************************************************
240    memstat=32
241
242    do indj=indmin,numbwf-1
243      if ((ldirect*wftime(indj).le.ldirect*itime).and. &
244           (ldirect*wftime(indj+1).gt.ldirect*itime)) then
245        memind(1)=1
246        if (lmpreader.or..not.lmp_use_reader) then
247          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
248            call readwind_ecmwf(indj,memind(1),uuh,vvh,wwh)
249          else
250            call readwind_gfs(indj,memind(1),uuh,vvh,wwh)
251          end if
252          call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn)
253          call calcpar(memind(1),uuh,vvh,pvh,metdata_format)
254          call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format)
255          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
256            call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh)
257          else
258            call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh)
259          end if
260          call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
261        end if
262        memtime(1)=wftime(indj)
263        memind(2)=2
264        if (lmpreader.or..not.lmp_use_reader) then
265          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
266            call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
267          else
268            call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh)
269          end if
270          call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
271          call calcpar(memind(2),uuh,vvh,pvh,metdata_format)
272          call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format)
273          if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
274            call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh)
275          else
276            call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh)
277          end if
278          call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
279        end if
280        memtime(2)=wftime(indj+1)
281! DEV: not used
282        if (.not.lmp_sync) memind(3)=3 ! indicate position of next read
283        nstop = 1
284        goto 60
285      endif
286    end do
28760  indmin=indj
288
289    mind3=memstat
290
291  endif
292
293  lwindinterv=abs(memtime(2)-memtime(1))
294
295  if (lwindinterv.gt.idiffmax) nstop=3
296
297end subroutine getfields
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG