source: flexpart.git/src/getfields.f90 @ 2eefa58

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file since 2eefa58 was 2eefa58, checked in by Espen Sollum ATMOS <eso@…>, 5 years ago

Added Ronas changes for inversion output
  • Property mode set to 100644
File size: 10.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,metdata_format)
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!   Unified ECMWF and GFS builds                                             *
45!   Marian Harustak, 12.5.2017                                               *
46!     - Added passing of metdata_format as it was needed by called routines  *
47!*****************************************************************************
48!                                                                            *
49! Variables:                                                                 *
50! lwindinterval [s]    time difference between the two wind fields read in   *
51! indj                 indicates the number of the wind field to be read in  *
52! indmin               remembers the number of wind fields already treated   *
53! memind(2)            pointer, on which place the wind fields are stored    *
54! memtime(2) [s]       times of the wind fields, which are kept in memory    *
55! itime [s]            current time since start date of trajectory calcu-    *
56!                      lation                                                *
57! nstop                > 0, if trajectory has to be terminated               *
58! nx,ny,nuvz,nwz       field dimensions in x,y and z direction               *
59! uu(0:nxmax,0:nymax,nuvzmax,2)   wind components in x-direction [m/s]       *
60! vv(0:nxmax,0:nymax,nuvzmax,2)   wind components in y-direction [m/s]       *
61! ww(0:nxmax,0:nymax,nwzmax,2)    wind components in z-direction [deltaeta/s]*
62! tt(0:nxmax,0:nymax,nuvzmax,2)   temperature [K]                            *
63! ps(0:nxmax,0:nymax,2)           surface pressure [Pa]                      *
64! metdata_format     format of metdata (ecmwf/gfs)                           *
65!                                                                            *
66! Constants:                                                                 *
67! idiffmax             maximum allowable time difference between 2 wind      *
68!                      fields                                                *
69!                                                                            *
70!*****************************************************************************
71
72  use par_mod
73  use com_mod
74  use class_gribfile
75
76  implicit none
77
78  integer :: indj,itime,nstop,memaux
79  integer :: metdata_format
80
81  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
82  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
83  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
84  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
85  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
86  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
87  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
88  real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
89! RLT added partial pressure water vapor
90  real :: pwater(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
91  integer :: kz, ix
92  character(len=100) :: rowfmt
93
94  integer :: indmin = 1
95
96
97! Check, if wind fields are available for the current time step
98!**************************************************************
99
100  nstop=0
101
102  if ((ldirect*wftime(1).gt.ldirect*itime).or. &
103       (ldirect*wftime(numbwf).lt.ldirect*itime)) then
104    write(*,*) 'FLEXPART WARNING: NO WIND FIELDS ARE AVAILABLE.'
105    write(*,*) 'A TRAJECTORY HAS TO BE TERMINATED.'
106    nstop=4
107    return
108  endif
109
110
111  if ((ldirect*memtime(1).le.ldirect*itime).and. &
112       (ldirect*memtime(2).gt.ldirect*itime)) then
113
114! The right wind fields are already in memory -> don't do anything
115!*****************************************************************
116
117    continue
118
119  else if ((ldirect*memtime(2).le.ldirect*itime).and. &
120       (memtime(2).ne.999999999)) then
121
122
123! Current time is after 2nd wind field
124! -> Resort wind field pointers, so that current time is between 1st and 2nd
125!***************************************************************************
126
127    memaux=memind(1)
128    memind(1)=memind(2)
129    memind(2)=memaux
130    memtime(1)=memtime(2)
131
132
133! Read a new wind field and store it on place memind(2)
134!******************************************************
135
136    do indj=indmin,numbwf-1
137      if (ldirect*wftime(indj+1).gt.ldirect*itime) then
138        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
139          call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
140        else
141          call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh)
142        end if
143        call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
144        call calcpar(memind(2),uuh,vvh,pvh,metdata_format)
145        call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format)
146        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
147          call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh)
148        else
149          call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh)
150        end if
151        call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
152        memtime(2)=wftime(indj+1)
153        nstop = 1
154        goto 40
155      endif
156    end do
15740  indmin=indj
158
159    if (WETBKDEP) then
160      call writeprecip(itime,memind(1))
161    endif
162
163  else
164
165! No wind fields, which can be used, are currently in memory
166! -> read both wind fields
167!***********************************************************
168
169    do indj=indmin,numbwf-1
170      if ((ldirect*wftime(indj).le.ldirect*itime).and. &
171           (ldirect*wftime(indj+1).gt.ldirect*itime)) then
172        memind(1)=1
173        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
174          call readwind_ecmwf(indj,memind(1),uuh,vvh,wwh)
175        else
176          call readwind_gfs(indj,memind(1),uuh,vvh,wwh)
177        end if
178        call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn)
179        call calcpar(memind(1),uuh,vvh,pvh,metdata_format)
180        call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format)
181        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
182          call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh)
183        else
184          call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh)
185        end if
186        call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
187        memtime(1)=wftime(indj)
188        memind(2)=2
189        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
190          call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
191        else
192          call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh)
193        end if
194        call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
195        call calcpar(memind(2),uuh,vvh,pvh,metdata_format)
196        call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format)
197        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
198          call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh)
199        else
200          call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh)
201        end if
202        call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
203        memtime(2)=wftime(indj+1)
204        nstop = 1
205        goto 60
206      endif
207    end do
20860  indmin=indj
209
210    if (WETBKDEP) then
211      call writeprecip(itime,memind(1))
212    endif
213
214  end if
215
216  ! RLT calculate dry air density
217  pwater=qv*prs/((r_air/r_water)*(1.-qv)+qv)
218  rho_dry=(prs-pwater)/(r_air*tt)
219
220  ! test density
221!  write(rowfmt,'(A,I6,A)') '(',nymax,'(E11.4,1X))'
222!  if(itime.eq.0) then
223!    open(500,file=path(2)(1:length(2))//'rho_dry.txt',status='replace',action='write')
224!    do kz=1,nzmax
225!      do ix=1,nxmax
226!        write(500,fmt=rowfmt) rho_dry(ix,:,kz,1)
227!      end do
228!    end do
229!    close(500)
230!    open(500,file=path(2)(1:length(2))//'rho.txt',status='replace',action='write')
231!    do kz=1,nzmax
232!      do ix=1,nxmax
233!        write(500,fmt=rowfmt) rho(ix,:,kz,1)
234!      end do
235!    end do
236!    close(500)
237!  endif
238
239  lwindinterv=abs(memtime(2)-memtime(1))
240
241  if (lwindinterv.gt.idiffmax) nstop=3
242
243end subroutine getfields
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG