source: branches/petra/src/getfields.f90 @ 37

Last change on this file since 37 was 37, checked in by pesei, 9 years ago

Wet dep quick fix and other small changes. Wet depo quick fix not final yet.
File size: 9.4 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)
23  !                       i     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  !  3/2015, PS: put in verbosity output
45  !
46  !*****************************************************************************
47  !                                                                            *
48  ! Variables:                                                                 *
49  ! lwindinterval [s]    time difference between the two wind fields read in   *
50  ! indj                 indicates the number of the wind field to be read in  *
51  ! indmin               remembers the number of wind fields already treated   *
52  ! memind(2)            pointer, on which place the wind fields are stored    *
53  ! memtime(2) [s]       times of the wind fields, which are kept in memory    *
54  ! itime [s]            current time since start date of trajectory calcu-    *
55  !                      lation                                                *
56  ! nstop                > 0, if trajectory has to be terminated               *
57  ! nx,ny,nuvz,nwz       field dimensions in x,y and z direction               *
58  ! uu(0:nxmax,0:nymax,nuvzmax,2)   wind components in x-direction [m/s]       *
59  ! vv(0:nxmax,0:nymax,nuvzmax,2)   wind components in y-direction [m/s]       *
60  ! ww(0:nxmax,0:nymax,nwzmax,2)    wind components in z-direction [deltaeta/s]*
61  ! tt(0:nxmax,0:nymax,nuvzmax,2)   temperature [K]                            *
62  ! ps(0:nxmax,0:nymax,2)           surface pressure [Pa]                      *
63  !                                                                            *
64  ! Constants:                                                                 *
65  ! idiffmax             maximum allowable time difference between 2 wind      *
66  !                      fields                                                *
67  !                                                                            *
68  !*****************************************************************************
69
70  use par_mod
71  use com_mod
72
73  implicit none
74
75  integer :: indj,itime,nstop,memaux
76
77  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
78  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
79  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
80  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
81  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
82  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
83  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
84  real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
85
86  integer :: indmin = 1
87
88  call verboutput(verbosity,1,'  getfields')
89
90  ! Check, if wind fields are available for the current time step
91  !**************************************************************
92
93  nstop=0
94
95  if ((ldirect*wftime(1).gt.ldirect*itime).or. &
96       (ldirect*wftime(numbwf).lt.ldirect*itime)) then
97    write(*,*) 'FLEXPART WARNING: NO WIND FIELDS ARE AVAILABLE.'
98    write(*,*) 'A TRAJECTORY HAS TO BE TERMINATED.'
99    nstop=4
100    return
101  endif
102
103
104  if ((ldirect*memtime(1).le.ldirect*itime).and. &
105       (ldirect*memtime(2).gt.ldirect*itime)) then
106
107  ! The right wind fields are already in memory -> don't do anything
108  !*****************************************************************
109
110    continue
111
112  else if ((ldirect*memtime(2).le.ldirect*itime).and. &
113       (memtime(2).ne.999999999)) then
114
115
116  ! Current time is after 2nd wind field
117  ! -> Resort wind field pointers, so that current time is between 1st and 2nd
118  !***************************************************************************
119
120    memaux=memind(1)
121    memind(1)=memind(2)
122    memind(2)=memaux
123    memtime(1)=memtime(2)
124
125
126  ! Read a new wind field and store it on place memind(2)
127  !******************************************************
128
129    do indj=indmin,numbwf-1
130       if (ldirect*wftime(indj+1).gt.ldirect*itime) then
131           call verboutput(verbosity,1,'    call readwind')
132          call readwind(indj+1,memind(2),uuh,vvh,wwh)
133           call verboutput(verbosity,1,'    call readwind_nests')
134          call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
135           call verboutput(verbosity,1,'    call calcpar')
136          call calcpar(memind(2),uuh,vvh,pvh)
137           call verboutput(verbosity,1,'    call calcpar_nests')
138          call calcpar_nests(memind(2),uuhn,vvhn,pvhn)
139           call verboutput(verbosity,1,'    call verttransform')
140          call verttransform(memind(2),uuh,vvh,wwh,pvh)
141           call verboutput(verbosity,1,'    call verttransform_nests')
142          call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
143          memtime(2)=wftime(indj+1)
144          nstop = 1
145          goto 40
146       endif
147    end do
148 40   indmin=indj
149
150  else
151
152  ! No wind fields, which can be used, are currently in memory
153  ! -> read both wind fields
154  !***********************************************************
155
156     do indj=indmin,numbwf-1
157        if ((ldirect*wftime(indj).le.ldirect*itime).and. &
158             (ldirect*wftime(indj+1).gt.ldirect*itime)) then
159           memind(1)=1
160            call verboutput(verbosity,1,'    call readwind')
161           call readwind(indj,memind(1),uuh,vvh,wwh)
162            call verboutput(verbosity,1,'    call readwind_nests')
163           call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn)
164            call verboutput(verbosity,1,'    call calcpar')
165           call calcpar(memind(1),uuh,vvh,pvh)
166            call verboutput(verbosity,1,'    call calcpar_nests')
167           call calcpar_nests(memind(1),uuhn,vvhn,pvhn)
168            call verboutput(verbosity,1,'    call verttransform')
169           call verttransform(memind(1),uuh,vvh,wwh,pvh)
170            call verboutput(verbosity,1,'    call verttransform_nests')
171           call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
172           memtime(1)=wftime(indj)
173           memind(2)=2
174            call verboutput(verbosity,1,'    call readwind')
175           call readwind(indj+1,memind(2),uuh,vvh,wwh)
176            call verboutput(verbosity,1,'    call readwind_nests')
177           call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
178            call verboutput(verbosity,1,'    call calcpar')
179           call calcpar(memind(2),uuh,vvh,pvh)
180            call verboutput(verbosity,1,'    call calcpar_nests')
181           call calcpar_nests(memind(2),uuhn,vvhn,pvhn)
182            call verboutput(verbosity,1,'    call verttransform')
183           call verttransform(memind(2),uuh,vvh,wwh,pvh)
184            call verboutput(verbosity,1,'    call verttransform_nests')
185           call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
186           memtime(2)=wftime(indj+1)
187           nstop = 1
188           goto 60
189        endif
190     end do
191 60   indmin=indj
192
193  endif
194
195  lwindinterv=abs(memtime(2)-memtime(1))
196
197  if (lwindinterv.gt.idiffmax) nstop=3
198
199  call verboutput(verbosity,1,'    leaving getfields')
200
201end subroutine getfields
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG