source: flexpart.git/src/getfields.f90

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

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