source: flexpart.git/src/getfields_with_dpdt.f90 @ 5f2e8f6

flexpart-noresm
Last change on this file since 5f2e8f6 was 5f2e8f6, checked in by Ignacio Pisso <Ignacio.Pisso@…>, 8 years ago

new flexpart noresm code in src
  • Property mode set to 100755
File size: 11.5 KB
Line 
1!**********************************************************************
2! Copyright 2016                                                      *
3! Andreas Stohl, Massimo Cassiani, Petra Seibert, A. Frank,           *
4! Gerhard Wotawa,  Caroline Forster, Sabine Eckhardt, John Burkhart,  *
5! Harald Sodemann                                                     *
6!                                                                     *
7! This file is part of FLEXPART-NorESM                                *
8!                                                                     *
9! FLEXPART-NorESM is free software: you can redistribute it           *
10! and/or modify                                                       *
11! it under the terms of the GNU General Public License as published by*
12! the Free Software Foundation, either version 3 of the License, or   *
13! (at your option) any later version.                                 *
14!                                                                     *
15! FLEXPART-NorESM is distributed in the hope that it will be useful,  *
16! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
17! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
18! GNU General Public License for more details.                        *
19!                                                                     *
20! You should have received a copy of the GNU General Public License   *
21! along with FLEXPART-NorESM.                                         *
22!  If not, see <http://www.gnu.org/licenses/>.                        *
23!**********************************************************************
24
25subroutine getfields(itime,nstop)
26  !                       i     o
27  !*****************************************************************************
28  !                                                                            *
29  !  This subroutine manages the 3 data fields to be kept in memory.           *
30  !  During the first time step of petterssen it has to be fulfilled that the  *
31  !  first data field must have |wftime|<itime, i.e. the absolute value of     *
32  !  wftime must be smaller than the absolute value of the current time in [s].*
33  !  The other 2 fields are the next in time after the first one.              *
34  !  Pointers (memind) are used, because otherwise one would have to resort the*
35  !  wind fields, which costs a lot of computing time. Here only the pointers  *
36  !  are resorted.                                                             *
37  !                                                                            *
38  !     Author: A. Stohl                                                       *
39  !     29 April 1994                                                          *
40  !                                                                            *
41  !*****************************************************************************
42  !  Changes, Bernd C. Krueger, Feb. 2001:
43  !   Variables tth,qvh,tthn,qvhn (on eta coordinates) in common block.
44  !   Function of nstop extended.
45  !*****************************************************************************
46  !  Modified by M. Cassiani    2016                                           *
47  !   - no nesting allowed in  FLEXPART-NorESM                                 *
48  !   - calcualate pressure derivative in time to be used in calcualting       *
49  !     vertical velocity with method_w=1                                      *
50  !   - two methods to calculate vertical velocity                             *
51  !                                                                            *
52  !*****************************************************************************
53  !                                                                            *
54  ! Variables:                                                                 *
55  ! lwindinterval [s]    time difference between the two wind fields read in   *
56  ! indj                 indicates the number of the wind field to be read in  *
57  ! indmin               remembers the number of wind fields already treated   *
58  ! memind(2)            pointer, on which place the wind fields are stored    *
59  ! memtime(2) [s]       times of the wind fields, which are kept in memory    *
60  ! itime [s]            current time since start date of trajectory calcu-    *
61  !                      lation                                                *
62  ! nstop                > 0, if trajectory has to be terminated               *
63  ! nx,ny,nuvz,nwz       field dimensions in x,y and z direction               *
64  ! uu(0:nxmax,0:nymax,nuvzmax,2)   wind components in x-direction [m/s]       *
65  ! vv(0:nxmax,0:nymax,nuvzmax,2)   wind components in y-direction [m/s]       *
66  ! ww(0:nxmax,0:nymax,nwzmax,2)    wind components in z-direction [deltaeta/s]*
67  ! tt(0:nxmax,0:nymax,nuvzmax,2)   temperature [K]                            *
68  ! ps(0:nxmax,0:nymax,2)           surface pressure [Pa]                      *
69  !                                                                            *
70  ! Constants:                                                                 *
71  ! idiffmax             maximum allowable time difference between 2 wind      *
72  !                      fields                                                *
73  !                                                                            *
74  !*****************************************************************************
75
76  use par_mod
77  use com_mod
78
79  implicit none
80
81  integer :: indj,itime,nstop,memaux
82
83  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
84  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
85  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
86  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
87  real :: uuhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
88  real :: vvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
89  real :: pvhn(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests)
90  real :: wwhn(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests)
91
92  integer :: indmin = 1
93
94  !*************************************************************
95  ! open(72,file='.\list_windfield.txt') !  for testing: mc
96  ! Check, if wind fields are available for the current time step
97  !**************************************************************
98
99  nstop=0
100
101  if ((ldirect*wftime(1).gt.ldirect*itime).or. &
102       (ldirect*wftime(numbwf).lt.ldirect*itime)) then
103    write(*,*) 'FLEXPART WARNING: NO WIND FIELDS ARE AVAILABLE.'
104    write(*,*) 'A TRAJECTORY HAS TO BE TERMINATED.'
105    nstop=4
106    return
107  endif
108
109
110  if ((ldirect*memtime(1).le.ldirect*itime).and. &
111       (ldirect*memtime(2).gt.ldirect*itime)) then
112
113  ! The right wind fields are already in memory -> don't do anything
114  !*****************************************************************
115
116    continue
117
118  else if ((ldirect*memtime(2).le.ldirect*itime).and. &
119       (memtime(2).ne.999999999)) then
120
121
122  ! Current time is after 2nd wind field
123  ! -> Resort wind field pointers, so that current time is between 1st and 2nd
124  !***************************************************************************
125
126    memaux=memind(1)
127    memind(1)=memind(2)
128    memind(2)=memaux
129    memtime(1)=memtime(2)
130
131
132  ! Read a new wind field and store it on place memind(2)
133  !******************************************************
134
135    do indj=indmin,numbwf-1
136       if (ldirect*wftime(indj+1).gt.ldirect*itime) then
137         
138          call read_delta_ps_intime(indj,1)   !load ps for previous (with respect to indj+1) time step
139          call read_delta_ps_intime(indj+2,2) !load ps for subsequent (with respect to indj+1) time step
140          call readwind(indj+1,memind(2),uuh,vvh,wwh)
141          !call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn) !nesting not active in NORESM version: comemnt by mc
142          call calcpar(memind(2),uuh,vvh,pvh)
143          !call calcpar_nests(memind(2),uuhn,vvhn,pvhn)         !nesting not active in NORESM version: comemnt by mc                   
144          if (method_w.eq.1) then !by mc 
145            call verttransform(memind(2),uuh,vvh,wwh,pvh,itime,indj+1)
146          else
147            call verttransform_omega(memind(2),uuh,vvh,wwh,pvh,itime,indj+1)
148          end if
149          !call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn) !nesting not active in NORESM version: comment by mc
150          memtime(2)=wftime(indj+1)
151          nstop = 1
152          goto 40
153       endif
154    end do
155 40   indmin=indj
156 
157  !*************************************************************************************************
158  !********************* for testing   by mc 
159  !  write(72,*)'getfields_with_dpdt: memind1    memtime1',memind(1), memtime(1) !for testin by mc
160  !  write(72,*)'getfields_with_dpdt: memind2    memtime2',memind(2), memtime(2) !for testin by mc
161  !**************************************************************************************************
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           call read_delta_ps_intime(indj-1,1) !load ps for previous (with respect to indj) time step
174           call read_delta_ps_intime(indj+1,2) !load ps for subsequent (with respect to indj) time step
175           call readwind(indj,memind(1),uuh,vvh,wwh)
176           !call readwind_nests(indj,memind(1),uuhn,vvhn,wwhn) !nesting not active in NORESM version: comment by mc
177           call calcpar(memind(1),uuh,vvh,pvh)
178           !call calcpar_nests(memind(1),uuhn,vvhn,pvhn) !nesting not active in NORESM version: comment by mc               
179           if (method_w.eq.1) then !added on 12-2-2016 by mc
180             call verttransform(memind(1),uuh,vvh,wwh,pvh,itime,indj)
181           else
182             call verttransform_omega(memind(1),uuh,vvh,wwh,pvh,itime,indj)   
183           end if
184           !call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn) !nesting not active in NORESM version: comment by mc
185           memtime(1)=wftime(indj)
186
187           memind(2)=2
188           call read_delta_ps_intime(indj,1)   !load ps for previous (with respect to indj+1) time step
189           call read_delta_ps_intime(indj+2,2) !load ps for subsequent (with respect to indj+1) time step
190           call readwind(indj+1,memind(2),uuh,vvh,wwh)
191           !call readwind_nests(indj+1,memind(2),uuhn,vvhn,wwhn) !nesting not active in NORESM version: comment by mc
192           call calcpar(memind(2),uuh,vvh,pvh)
193           !call calcpar_nests(memind(2),uuhn,vvhn,pvhn) !nesting not active in NORESM version: commemnt by mc
194           if (method_w.eq.1) then !added on 12-2-2016 by mc
195             call verttransform(memind(2),uuh,vvh,wwh,pvh,itime,indj+1)
196           else
197             call verttransform_omega(memind(2),uuh,vvh,wwh,pvh,itime,indj+1)   
198           end if         
199           !call verttransform_nests(memind(2),uuhn,vvhn,wwhn,pvhn) !nesting not active in NORESM version: comment by mc
200           memtime(2)=wftime(indj+1)
201           nstop = 1
202           goto 60
203        endif
204     end do
205 60   indmin=indj
206
207!******************************************************************************************************
208!********************* for testing  by mc   
209!     write(72,*)'getfields_with_dpdt: memind1    memtime1',memind(1), memtime(1) !for testin by mc
210!     write(72,*)'getfields_with_dpdt: memind2    memtime2',memind(2), memtime(2) !for testin by mc
211!*******************************************************************************************************
212  endif
213
214  lwindinterv=abs(memtime(2)-memtime(1))
215
216  if (lwindinterv.gt.idiffmax) nstop=3
217
218end subroutine getfields
Note: See TracBrowser for help on using the repository browser.
hosted by ZAMG