source: flexpart.git/src/getfields.f90 @ 3481cc1

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bug
Last change on this file since 3481cc1 was 3481cc1, checked in by Ignacio Pisso <ip@…>, 4 years ago

move license from headers to a different file
  • Property mode set to 100644
File size: 9.0 KB
Line 
1subroutine getfields(itime,nstop,metdata_format)
2!                       i     o
3!*****************************************************************************
4!                                                                            *
5!  This subroutine manages the 3 data fields to be kept in memory.           *
6!  During the first time step of petterssen it has to be fulfilled that the  *
7!  first data field must have |wftime|<itime, i.e. the absolute value of     *
8!  wftime must be smaller than the absolute value of the current time in [s].*
9!  The other 2 fields are the next in time after the first one.              *
10!  Pointers (memind) are used, because otherwise one would have to resort the*
11!  wind fields, which costs a lot of computing time. Here only the pointers  *
12!  are resorted.                                                             *
13!                                                                            *
14!     Author: A. Stohl                                                       *
15!                                                                            *
16!     29 April 1994                                                          *
17!                                                                            *
18!*****************************************************************************
19!  Changes, Bernd C. Krueger, Feb. 2001:                                     *
20!   Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block.        *
21!   Function of nstop extended.                                              *
22!                                                                            *
23!   Unified ECMWF and GFS builds                                             *
24!   Marian Harustak, 12.5.2017                                               *
25!     - Added passing of metdata_format as it was needed by called routines  *
26!*****************************************************************************
27!                                                                            *
28! Variables:                                                                 *
29! lwindinterval [s]    time difference between the two wind fields read in   *
30! indj                 indicates the number of the wind field to be read in  *
31! indmin               remembers the number of wind fields already treated   *
32! memind(2)            pointer, on which place the wind fields are stored    *
33! memtime(2) [s]       times of the wind fields, which are kept in memory    *
34! itime [s]            current time since start date of trajectory calcu-    *
35!                      lation                                                *
36! nstop                > 0, if trajectory has to be terminated               *
37! nx,ny,nuvz,nwz       field dimensions in x,y and z direction               *
38! uu(0:nxmax,0:nymax,nuvzmax,2)   wind components in x-direction [m/s]       *
39! vv(0:nxmax,0:nymax,nuvzmax,2)   wind components in y-direction [m/s]       *
40! ww(0:nxmax,0:nymax,nwzmax,2)    wind components in z-direction [deltaeta/s]*
41! tt(0:nxmax,0:nymax,nuvzmax,2)   temperature [K]                            *
42! ps(0:nxmax,0:nymax,2)           surface pressure [Pa]                      *
43! metdata_format     format of metdata (ecmwf/gfs)                           *
44!                                                                            *
45! Constants:                                                                 *
46! idiffmax             maximum allowable time difference between 2 wind      *
47!                      fields                                                *
48!                                                                            *
49!*****************************************************************************
50
51  use par_mod
52  use com_mod
53  use class_gribfile
54
55  implicit none
56
57  integer :: indj,itime,nstop,memaux
58  integer :: metdata_format
59
60  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
61  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
62  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
63  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
64  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
65  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
66  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
67  real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
68! RLT added partial pressure water vapor
69  real :: pwater(0:nxmax-1,0:nymax-1,nzmax,numwfmem)
70  integer :: kz, ix
71  character(len=100) :: rowfmt
72
73  integer :: indmin = 1
74
75
76! Check, if wind fields are available for the current time step
77!**************************************************************
78
79  nstop=0
80
81  if ((ldirect*wftime(1).gt.ldirect*itime).or. &
82       (ldirect*wftime(numbwf).lt.ldirect*itime)) then
83    write(*,*) 'FLEXPART WARNING: NO WIND FIELDS ARE AVAILABLE.'
84    write(*,*) 'A TRAJECTORY HAS TO BE TERMINATED.'
85    nstop=4
86    return
87  endif
88
89
90  if ((ldirect*memtime(1).le.ldirect*itime).and. &
91       (ldirect*memtime(2).gt.ldirect*itime)) then
92
93! The right wind fields are already in memory -> don't do anything
94!*****************************************************************
95
96    continue
97
98  else if ((ldirect*memtime(2).le.ldirect*itime).and. &
99       (memtime(2).ne.999999999)) then
100
101
102! Current time is after 2nd wind field
103! -> Resort wind field pointers, so that current time is between 1st and 2nd
104!***************************************************************************
105
106    memaux=memind(1)
107    memind(1)=memind(2)
108    memind(2)=memaux
109    memtime(1)=memtime(2)
110
111
112! Read a new wind field and store it on place memind(2)
113!******************************************************
114
115    do indj=indmin,numbwf-1
116      if (ldirect*wftime(indj+1).gt.ldirect*itime) then
117        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
118          call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
119        else
120          call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh)
121        end if
122        call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
123        call calcpar(memind(2),uuh,vvh,pvh,metdata_format)
124        call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format)
125        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
126          call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh)
127        else
128          call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh)
129        end if
130        call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
131        memtime(2)=wftime(indj+1)
132        nstop = 1
133        goto 40
134      endif
135    end do
13640  indmin=indj
137
138    if (WETBKDEP) then
139      call writeprecip(itime,memind(1))
140    endif
141
142  else
143
144! No wind fields, which can be used, are currently in memory
145! -> read both wind fields
146!***********************************************************
147
148    do indj=indmin,numbwf-1
149      if ((ldirect*wftime(indj).le.ldirect*itime).and. &
150           (ldirect*wftime(indj+1).gt.ldirect*itime)) then
151        memind(1)=1
152        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
153          call readwind_ecmwf(indj,memind(1),uuh,vvh,wwh)
154        else
155          call readwind_gfs(indj,memind(1),uuh,vvh,wwh)
156        end if
157        call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn)
158        call calcpar(memind(1),uuh,vvh,pvh,metdata_format)
159        call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format)
160        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
161          call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh)
162        else
163          call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh)
164        end if
165        call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
166        memtime(1)=wftime(indj)
167        memind(2)=2
168        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
169          call readwind_ecmwf(indj+1,memind(2),uuh,vvh,wwh)
170        else
171          call readwind_gfs(indj+1,memind(2),uuh,vvh,wwh)
172        end if
173        call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn)
174        call calcpar(memind(2),uuh,vvh,pvh,metdata_format)
175        call calcpar_nests(memind(2),uuhn,vvhn,pvhn,metdata_format)
176        if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
177          call verttransform_ecmwf(memind(2),uuh,vvh,wwh,pvh)
178        else
179          call verttransform_gfs(memind(2),uuh,vvh,wwh,pvh)
180        end if
181        call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn)
182        memtime(2)=wftime(indj+1)
183        nstop = 1
184        goto 60
185      endif
186    end do
18760  indmin=indj
188
189    if (WETBKDEP) then
190      call writeprecip(itime,memind(1))
191    endif
192
193  end if
194
195  ! RLT calculate dry air density
196  pwater=qv*prs/((r_air/r_water)*(1.-qv)+qv)
197  rho_dry=(prs-pwater)/(r_air*tt)
198
199  ! test density
200!  write(rowfmt,'(A,I6,A)') '(',nymax,'(E11.4,1X))'
201!  if(itime.eq.0) then
202!    open(500,file=path(2)(1:length(2))//'rho_dry.txt',status='replace',action='write')
203!    do kz=1,nzmax
204!      do ix=1,nxmax
205!        write(500,fmt=rowfmt) rho_dry(ix,:,kz,1)
206!      end do
207!    end do
208!    close(500)
209!    open(500,file=path(2)(1:length(2))//'rho.txt',status='replace',action='write')
210!    do kz=1,nzmax
211!      do ix=1,nxmax
212!        write(500,fmt=rowfmt) rho(ix,:,kz,1)
213!      end do
214!    end do
215!    close(500)
216!  endif
217
218  lwindinterv=abs(memtime(2)-memtime(1))
219
220  if (lwindinterv.gt.idiffmax) nstop=3
221
222end subroutine getfields
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG