1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
|
subroutine readcet(error)
C o
***********************************************************************
* *
* TRAJECTORY MODEL SUBROUTINE READCET *
* *
***********************************************************************
* *
* AUTHOR: A. STOHL *
* DATE: 1999-01-11 *
* *
* Update: Vertical coordinate can now be given in meters above sea *
* level, meters above ground, and in hPa. *
* *
***********************************************************************
* *
* DESCRIPTION: *
* *
* READING OF TRAJECTORY STARTING/ENDING POINTS FROM DATA FILE *
* *
* LINE a line of text *
* NUMPOINT number of trajectory starting/ending points *
* XPOINT(maxpoint) x-coordinates of starting/ending points *
* YPOINT(maxpoint) y-coordinates of starting/ending points *
* ZPOINT(maxpoint) z-coordinates of starting/ending points *
* KINDZ(maxpoint) kind of z coordinate (1:masl, 2:magl, 3:hPa) *
* KIND(maxpoint) kind of trajectory *
* COMPOINT(maxpoint) comment for trajectory output *
* *
***********************************************************************
*
include 'includepar'
include 'includecom'
integer isum
logical error,old
integer kindhelp,kindzhelp
real xhelpleft,xhelpright,yhelpleft,yhelpright,deltax
real deltay,zhelplower,zhelpupper,deltaz,xi,yj,zk
character*45 comhelp,line
error=.false.
*
* OPENING OF FILE 'STARTCET'
*
open(unitpoin,file=path(1)(1:len(1))//'STARTCET',
+status='old',err=999)
C Check the format of the STARTCET file (either in free format,
C or using formatted mask)
C Use of formatted mask is assumed if line 28 contains the word '____'
**********************************************************************
call skplin(27,unitpoin)
read (unitpoin,901) line
901 format (a)
if (index(line,'____') .eq. 0) then
old = .false.
else
old = .true.
endif
rewind(unitpoin)
*
* READING OF STARTING/ENDING POINTS OF TRAJECTORIES
*
call skplin(26,unitpoin)
read(unitpoin,*,err=998,end=998) xhelpleft
if (old) call skplin(2,unitpoin)
read(unitpoin,*,err=998,end=998) yhelpleft
if (old) call skplin(2,unitpoin)
read(unitpoin,*,err=998,end=998) xhelpright
if (old) call skplin(2,unitpoin)
read(unitpoin,*,err=998,end=998) yhelpright
if (old) call skplin(2,unitpoin)
read(unitpoin,*,err=998,end=998) deltax
if (old) call skplin(2,unitpoin)
read(unitpoin,*,err=998,end=998) deltay
if (old) call skplin(2,unitpoin)
read(unitpoin,*,err=998,end=998) kindhelp
if (old) call skplin(2,unitpoin)
read(unitpoin,*,err=998,end=998) kindzhelp
if (old) call skplin(2,unitpoin)
read(unitpoin,*,err=998,end=998) zhelplower
if (old) call skplin(2,unitpoin)
read(unitpoin,*,err=998,end=998) zhelpupper
if (old) call skplin(2,unitpoin)
read(unitpoin,*,err=998,end=998) deltaz
if (old) call skplin(2,unitpoin)
if (old) then
read(unitpoin,'(a40)',err=998,end=998) comhelp(1:40)
else
read(unitpoin,*,err=998,end=998) comhelp(1:40)
endif
if (old) call skplin(1,unitpoin)
close(unitpoin)
1 format(1x,a17,i6,a28)
if (zhelpupper.lt.zhelplower) stop 'Upper level must be greater
+than lower level'
C Forbid mixing layer trajectories
**********************************
if (kindhelp.eq.3) then
write(*,*) '### Mixing layer trajectories not allowed ###'
write(*,*) '### for CET calculations. Select a ###'
write(*,*) '### different trajectory type! ###'
error=.true.
return
endif
C Check field dimension
***********************
isum=int(((xhelpright-xhelpleft)/deltax+1.)*((yhelpright-
+yhelpleft)/deltay+1.)*(abs(zhelpupper-zhelplower)/deltaz +1.))
if (isum.gt.maxtra) then
write(*,*) ' #### TRAJECTORY MODEL ERROR! MODEL CAN NOT #### '
write(*,*) ' #### KEEP SO MANY TRAJECTORIES IN MEMORY. #### '
write(*,1) ' #### CURRENTLY: ',isum,' ####
+ '
write(*,1) ' #### MAXIMUM: ',maxtra,' ####
+ '
write(*,*) ' #### REDUCE CET DOMAIN OR RESOLUTION. #### '
error=.true.
endif
C Initialize CET starting points
********************************
isum=0
do 10 xi=xhelpleft,xhelpright,deltax
do 10 yj=yhelpleft,yhelpright,deltay
do 10 zk=zhelplower,zhelpupper,deltaz
isum=isum+1
xpoint(isum)=xi
ypoint(isum)=yj
kind(isum)=kindhelp
kindz(isum)=kindzhelp
zpoint(isum)=zk
compoint(isum)(1:40)=comhelp(1:40)
C Conversion from hPa to Pa, if z coordinate is given in pressure units
if (kindz(isum).eq.3) zpoint(isum)=zpoint(isum)*100.
10 continue
numpoint=isum
return
998 error=.true.
write(*,*) '#### TRAJECTORY MODEL SUBROUTINE READCET: '//
& '#### '
write(*,*) '#### FATAL ERROR - FILE "STARTCET" IS '//
& '#### '
write(*,*) '#### CORRUPT. PLEASE CHECK YOUR INPUTS FOR '//
& '#### '
write(*,*) '#### MISTAKES OR GET A NEW "STARTPOINTS"- '//
& '#### '
write(*,*) '#### FILE ... '//
& '#### '
return
999 error=.true.
write(*,*)
write(*,*) ' ###########################################'//
& '###### '
write(*,*) ' TRAJECTORY MODEL SUBROUTINE READCET:'
write(*,*)
write(*,*) ' FATAL ERROR - FILE STARTCET IS NOT AVAILABLE'
write(*,*) ' OR YOU ARE NOT PERMITTED FOR ANY ACCESS'
write(*,*) ' ###########################################'//
& '###### '
return
end
|
