Subversion Repositories lagranto.ecmwf

Rev

Rev 3 | Only display areas with differences | Ignore whitespace | Details | Blame | Last modification | View Log | RSS feed

Rev 3 Rev 44
1 
c     -------------------------------------------------------------------------
 1 
c     -------------------------------------------------------------------------
2 
c     Potential temperature (TH)
 2 
c     Potential temperature (TH)
3 
c     -------------------------------------------------------------------------
 3 
c     -------------------------------------------------------------------------
4 
 
 4 
 
5 
      subroutine calc_TH (pt,t,p)
 5 
      subroutine calc_TH (pt,t,p)
6
6
7 
      implicit none
 7 
      implicit none
8
8
9 
c     Argument declaration
 9 
c     Argument declaration
10 
      real  pt     ! Potential temperature [K]
 10 
      real  pt     ! Potential temperature [K]
11 
      real  t      ! Temperature [either in C or in K]
 11 
      real  t      ! Temperature [either in C or in K]
12 
      real  p      ! Pressure [hPa]
 12 
      real  p      ! Pressure [hPa]
13 
 
 13 
 
14 
c     Physical parameters
 14 
c     Physical parameters
15 
      real      rdcp,tzero
 15 
      real      rdcp,tzero
16 
      data      rdcp,tzero /0.286,273.15/
 16 
      data      rdcp,tzero /0.286,273.15/
17
17
18 
c     Calculation - distinction between temperature in C or in K
 18 
c     Calculation - distinction between temperature in C or in K
19 
      if (t.lt.100.) then
 19 
      if (t.lt.100.) then
20 
         pt = (t+tzero) * ( (1000./p)**rdcp )
20 
         pt = (t+tzero) * ( (1000./p)**rdcp )
21 
      else
 21 
      else
22 
         pt = t * ( (1000./p)**rdcp )
 22 
         pt = t * ( (1000./p)**rdcp )
23 
      endif
 23 
      endif
24 
 
 24 
 
25 
      end
 25 
      end
26 
 
 26 
 
27 
c     -------------------------------------------------------------------------
 27 
c     -------------------------------------------------------------------------
28 
c     Density (RHO)
 28 
c     Density (RHO)
29 
c     -------------------------------------------------------------------------
 29 
c     -------------------------------------------------------------------------
30 
 
 30 
 
31 
      subroutine calc_RHO (rho,t,p)
 31 
      subroutine calc_RHO (rho,t,p)
32
32
33 
      implicit none
 33 
      implicit none
34 
 
 34 
 
35 
c     Argument declaration
 35 
c     Argument declaration
36 
      real  rho    ! Density [kg/m^3]
 36 
      real  rho    ! Density [kg/m^3]
37 
      real  t      ! Temperature [either in C or in K]
 37 
      real  t      ! Temperature [either in C or in K]
38 
      real  p      ! Pressure [hPa]
 38 
      real  p      ! Pressure [hPa]
39 
 
 39 
 
40 
c     Physical parameters
 40 
c     Physical parameters
41 
      real      rd,tzero
 41 
      real      rd,tzero
42 
      data      rd,tzero /287.05,273.15/
 42 
      data      rd,tzero /287.05,273.15/
43 
 
 43 
 
44 
c     Auxiliary variables
 44 
c     Auxiliary variables
45 
      real      tk
 45 
      real      tk
46
46
47 
c     Calculation - distinction between temperature in C or in K
 47 
c     Calculation - distinction between temperature in C or in K
48 
      if (t.lt.100.) then
 48 
      if (t.lt.100.) then
49 
         tk = t + tzero
 49 
         tk = t + tzero
50 
      else
 50 
      else
51 
         tk = t
 51 
         tk = t
52 
      endif
 52 
      endif
53 
 
 53 
 
54 
      rho = 100.*p/( tk * rd )
54 
      rho = 100.*p/( tk * rd )
55 
 
 55 
 
56 
      end
 56 
      end
57 
 
 57 
 
58 
c     -------------------------------------------------------------------------
 58 
c     -------------------------------------------------------------------------
59 
c     Relative humidity (RH)
 59 
c     Relative humidity (RH)
60 
c     -------------------------------------------------------------------------
 60 
c     -------------------------------------------------------------------------
61 
 
 61 
 
62 
      subroutine calc_RH (rh,t,p,q)
 62 
      subroutine calc_RH (rh,t,p,q)
63
63
64 
      implicit none
 64 
      implicit none
65 
 
 65 
 
66 
c     Argument declaration
 66 
c     Argument declaration
67 
      real  rh     ! Relative humidity [%]
 67 
      real  rh     ! Relative humidity [%]
68 
      real  t      ! Temperature [either in C or in K]
 68 
      real  t      ! Temperature [either in C or in K]
69 
      real  p      ! Pressure [hPa]
 69 
      real  p      ! Pressure [hPa]
70 
      real  q      ! Specific humidity [kg/kg]
 70 
      real  q      ! Specific humidity [kg/kg]
71 
 
 71 
 
72 
c     Physical parameters
 72 
c     Physical parameters
73 
      real      rdcp,tzero
 73 
      real      rdcp,tzero
74 
      data      rdcp,tzero /0.286,273.15/
 74 
      data      rdcp,tzero /0.286,273.15/
75 
      real      b1,b2w,b3,b4w,r,rd
 75 
      real      b1,b2w,b3,b4w,r,rd
76 
      data      b1,b2w,b3,b4w,r,rd /6.1078, 17.2693882, 273.16, 35.86,
 76 
      data      b1,b2w,b3,b4w,r,rd /6.1078, 17.2693882, 273.16, 35.86,
77 
     &                              287.05, 461.51/
 77 
     &                              287.05, 461.51/
78 
 
 78 
 
79 
c     Auxiliary variables
 79 
c     Auxiliary variables
80 
      real      ge
 80 
      real      ge
81 
      real      gqd
 81 
      real      gqd
82 
      real      tc
 82 
      real      tc
83 
      real      pp,qk
 83 
      real      pp,qk
84 
 
 84 
 
85 
c     Calculation - distinction between temperature in C or in K
 85 
c     Calculation - distinction between temperature in C or in K
86 
      if (t.gt.100.) then
 86 
      if (t.gt.100.) then
87 
         tc = t - tzero
 87 
         tc = t - tzero
88 
      else
 88 
      else
89 
         tc = t
 89 
         tc = t
90 
      endif
 90 
      endif
91 
      qk = q
 91 
      qk = q
92 
 
 92 
 
93 
      ge  = b1*exp(b2w*tc/(tc+b3-b4w))
 93 
      ge  = b1*exp(b2w*tc/(tc+b3-b4w))
94 
      gqd = r/rd*ge/(p-(1.-r/rd)*ge)
 94 
      gqd = r/rd*ge/(p-(1.-r/rd)*ge)
95 
      rh  = 100.*qk/gqd
 95 
      rh  = 100.*qk/gqd
96 
 
 96 
 
97 
      end
 97 
      end
98 
 
 98 
 
99 
c     -------------------------------------------------------------------------
 99 
c     -------------------------------------------------------------------------
100 
c     Equivalent potential temperature (THE)
 100 
c     Equivalent potential temperature (THE)
101 
c     -------------------------------------------------------------------------
 101 
c     -------------------------------------------------------------------------
102 
 
 102 
 
103 
      subroutine calc_THE (the,t,p,q)
 103 
      subroutine calc_THE (the,t,p,q)
104 
 
 104 
 
105 
      implicit none
 105 
      implicit none
106 
 
 106 
 
107 
c     Argument declaration
 107 
c     Argument declaration
108 
      real  the    ! Equivalent potential temperature [K]
 108 
      real  the    ! Equivalent potential temperature [K]
109 
      real  t      ! Temperature [either in C or in K]
 109 
      real  t      ! Temperature [either in C or in K]
110 
      real  p      ! Pressure [hPa]
 110 
      real  p      ! Pressure [hPa]
111 
      real  q      ! Specific humidity [kg/kg]
 111 
      real  q      ! Specific humidity [kg/kg]
112 
 
 112 
 
113 
c     Physical parameters
113 
c     Physical parameters
114 
      real     rdcp,tzero
 114 
      real     rdcp,tzero
115 
      data     rdcp,tzero /0.286,273.15/
 115 
      data     rdcp,tzero /0.286,273.15/
116 
 
 116 
 
117 
c     Auxiliary variables
117 
c     Auxiliary variables
118 
      real     tk,qk
 118 
      real     tk,qk
119
119
120 
c     Calculation - distinction between temperature in C or in K
 120 
c     Calculation - distinction between temperature in C or in K
121 
      if (t.lt.100.) then
 121 
      if (t.lt.100.) then
122 
         tk = t + tzero
 122 
         tk = t + tzero
123 
      else
 123 
      else
124 
         tk = t
 124 
         tk = t
125 
      endif
 125 
      endif
126 
      qk = q
 126 
      qk = q
127
127
128 
      the = tk*(1000./p)
 128 
      the = tk*(1000./p)
129 
     +      **(0.2854*(1.0-0.28*qk))*exp(
 129 
     +      **(0.2854*(1.0-0.28*qk))*exp(
130 
     +      (3.376/(2840.0/(3.5*alog(tk)-alog(
 130 
     +      (3.376/(2840.0/(3.5*alog(tk)-alog(
131 
     +      100.*p*max(1.0E-10,qk)/(0.622+0.378*
 131 
     +      100.*p*max(1.0E-10,qk)/(0.622+0.378*
132 
     +      q))-0.1998)+55.0)-0.00254)*1.0E3*
 132 
     +      q))-0.1998)+55.0)-0.00254)*1.0E3*
133 
     +      max(1.0E-10,qk)*(1.0+0.81*qk))
 133 
     +      max(1.0E-10,qk)*(1.0+0.81*qk))
134 
 
 134 
 
135 
      end
 135 
      end
136 
 
 136 
 
137 
c     -------------------------------------------------------------------------
 137 
c     -------------------------------------------------------------------------
138 
c     Latent heating rate (LHR)
 138 
c     Latent heating rate (LHR)
139 
c     -------------------------------------------------------------------------
 139 
c     -------------------------------------------------------------------------
140 
 
 140 
 
141 
      subroutine calc_LHR (lhr,t,p,q,omega,rh)
 141 
      subroutine calc_LHR (lhr,t,p,q,omega,rh)
142 
 
 142 
 
143 
      implicit none
143 
      implicit none
144 
 
 144 
 
145 
c     Argument declaration
 145 
c     Argument declaration
146 
      real  lhr    ! Latent heating rate [K/6h]
 146 
      real  lhr    ! Latent heating rate [K/6h]
147 
      real  t      ! Temperature [either in C or in K]
 147 
      real  t      ! Temperature [either in C or in K]
148 
      real  p      ! Pressure [hPa]
 148 
      real  p      ! Pressure [hPa]
149 
      real  q      ! Specific humidity [kg/kg]
 149 
      real  q      ! Specific humidity [kg/kg]
150 
      real  omega  ! Vertical velocity [Pa/s]
 150 
      real  omega  ! Vertical velocity [Pa/s]
151 
      real  rh     ! Relative humidity [%]
151 
      real  rh     ! Relative humidity [%]
152 
 
 152 
 
153 
c     Physical parameters
153 
c     Physical parameters
154 
      real      p0,kappa,tzero
 154 
      real      p0,kappa,tzero
155 
      data      p0,kappa,tzero /1000.,0.286,273.15/
 155 
      data      p0,kappa,tzero /1000.,0.286,273.15/
156 
      real      blog10,cp,r,lw,eps
 156 
      real      blog10,cp,r,lw,eps
157 
      data      blog10,cp,r,lw,eps /.08006,1004.,287.,2.5e+6,0.622/
 157 
      data      blog10,cp,r,lw,eps /.08006,1004.,287.,2.5e+6,0.622/
158 
 
 158 
 
159 
c     Auxiliary variables
 159 
c     Auxiliary variables
160 
      real      tk
 160 
      real      tk
161 
      real      qk
 161 
      real      qk
162 
      real      tt
 162 
      real      tt
163 
      real      esat,c
 163 
      real      esat,c
164 
 
 164 
 
165 
c     Calculation - distinction between temperature in C or in K
 165 
c     Calculation - distinction between temperature in C or in K
166 
      if (t.lt.100.) then
 166 
      if (t.lt.100.) then
167 
         tk = t + tzero
 167 
         tk = t + tzero
168 
      else
 168 
      else
169 
         tk = t
 169 
         tk = t
170 
      endif
 170 
      endif
171 
      qk = q
 171 
      qk = q
172 
 
 172 
 
173 
      if (rh.lt.80.) then
173 
      if (rh.lt.80.) then
174 
         lhr = 0.
174 
         lhr = 0.
175 
      else if (omega.gt.0.) then
175 
      else if (omega.gt.0.) then
176 
         lhr = 0.
 176 
         lhr = 0.
177 
      else
 177 
      else
178 
         c   = lw/cp*eps*blog10*esat(tk)/p
 178 
         c   = lw/cp*eps*blog10*esat(tk)/p
179 
         tt  = (tk*(p0/p)**kappa)
179 
         tt  = (tk*(p0/p)**kappa)
180 
         lhr = 21600.*
180 
         lhr = 21600.*
181 
     >        (1.-exp(.2*(80.-rh)))
181 
     >        (1.-exp(.2*(80.-rh)))
182 
     >        *(-c*kappa*tt*omega/(100.*p))/(1.+c)
182 
     >        *(-c*kappa*tt*omega/(100.*p))/(1.+c)
183
183
184 
      endif
 184 
      endif
185
185
186 
      end
 186 
      end
187 
 
 187 
 
188 
c     -------------------------------------------------------------------------
 188 
c     -------------------------------------------------------------------------
189 
c     Wind speed (VEL)
 189 
c     Wind speed (VEL)
190 
c     -------------------------------------------------------------------------
 190 
c     -------------------------------------------------------------------------
191 
 
 191 
 
192 
      subroutine calc_VEL (vel,u,v)
 192 
      subroutine calc_VEL (vel,u,v)
193 
 
 193 
 
194 
      implicit none
 194 
      implicit none
195
195
196 
c     Argument declaration
 196 
c     Argument declaration
197 
      real  vel    ! Wind speed [m/s]
 197 
      real  vel    ! Wind speed [m/s]
198 
      real  u      ! Zonal wind component [m/s]
 198 
      real  u      ! Zonal wind component [m/s]
199 
      real  v      ! Meridional wind component [m/s]
 199 
      real  v      ! Meridional wind component [m/s]
200 
 
 200 
 
201 
      vel = sqrt ( u*u + v*v )
 201 
      vel = sqrt ( u*u + v*v )
202 
 
 202 
 
203 
      end
 203 
      end
204 
 
 204 
 
205 
c     -------------------------------------------------------------------------
 205 
c     -------------------------------------------------------------------------
206 
c     Wind direction (DIR)
 206 
c     Wind direction (DIR)
207 
c     -------------------------------------------------------------------------
 207 
c     -------------------------------------------------------------------------
208 
 
 208 
 
209 
      subroutine calc_DIR (dir,u,v)
 209 
      subroutine calc_DIR (dir,u,v)
210 
 
 210 
 
211 
      implicit none
 211 
      implicit none
212
212
213 
c     Argument declaration
 213 
c     Argument declaration
214 
      real  dir    ! Wind direction [deg]
 214 
      real  dir    ! Wind direction [deg]
215 
      real  u      ! Zonal wind component [m/s]
 215 
      real  u      ! Zonal wind component [m/s]
216 
      real  v      ! Meridional wind component [m/s]
 216 
      real  v      ! Meridional wind component [m/s]
217 
 
 217 
 
218 
      call getangle(1.,0.,u,v,dir)
 218 
      call getangle(1.,0.,u,v,dir)
219 
 
 219 
 
220 
      end
 220 
      end
221 
 
 221 
 
222 
c     -------------------------------------------------------------------------
 222 
c     -------------------------------------------------------------------------
223 
c     Zonal derivative of U (DUDX)
 223 
c     Zonal derivative of U (DUDX)
224 
c     -------------------------------------------------------------------------
 224 
c     -------------------------------------------------------------------------
225 
 
 225 
 
226 
      subroutine calc_DUDX (dudx,u1,u0,lat)
 226 
      subroutine calc_DUDX (dudx,u1,u0,lat)
227 
 
 227 
 
228 
      implicit none
 228 
      implicit none
229
229
230 
c     Argument declaration
 230 
c     Argument declaration
231 
      real  dudx   ! Derivative of U in zonal direction [s^-1]
 231 
      real  dudx   ! Derivative of U in zonal direction [s^-1]
232 
      real  u1     ! U @ LON + 1 DLON [m/s]
 232 
      real  u1     ! U @ LON + 1 DLON [m/s]
233 
      real  u0     ! U @ LON - 1 DLON [m/s]
 233 
      real  u0     ! U @ LON - 1 DLON [m/s]
234 
      real  lat    ! Latitude [deg]
 234 
      real  lat    ! Latitude [deg]
235 
 
 235 
 
236 
c     Physical parameters
 236 
c     Physical parameters
237 
      real      pi180
237 
      real      pi180
238 
      parameter (pi180=31.14159/180.)
 238 
      parameter (pi180=31.14159/180.)
239 
      real      deltay
 239 
      real      deltay
240 
      parameter (deltay =1.11e5)
 240 
      parameter (deltay =1.11e5)
241
241
242 
      dudx = (u1-u0) / ( 2. * deltay * cos(pi180 * lat) )
 242 
      dudx = (u1-u0) / ( 2. * deltay * cos(pi180 * lat) )
243 
 
 243 
 
244 
      end
 244 
      end
245 
 
 245 
 
246 
c     -------------------------------------------------------------------------
 246 
c     -------------------------------------------------------------------------
247 
c     Zonal derivative of V (DVDX)
 247 
c     Zonal derivative of V (DVDX)
248 
c     -------------------------------------------------------------------------
 248 
c     -------------------------------------------------------------------------
249 
 
 249 
 
250 
      subroutine calc_DVDX (dvdx,v1,v0,lat)
 250 
      subroutine calc_DVDX (dvdx,v1,v0,lat)
251
251
252 
c     Argument declaration
 252 
c     Argument declaration
253 
      real  dvdx   ! Derivative of V in zonal direction [s^-1]
 253 
      real  dvdx   ! Derivative of V in zonal direction [s^-1]
254 
      real  v1     ! V @ LON + 1 DLON [m/s]
 254 
      real  v1     ! V @ LON + 1 DLON [m/s]
255 
      real  v0     ! V @ LON - 1 DLON [m/s]
 255 
      real  v0     ! V @ LON - 1 DLON [m/s]
256 
      real  lat    ! Latitude [deg]
 256 
      real  lat    ! Latitude [deg]
257 
 
 257 
 
258 
c     Physical parameters
 258 
c     Physical parameters
259 
      real      pi180
259 
      real      pi180
260 
      parameter (pi180=3.14159/180.)
 260 
      parameter (pi180=3.14159/180.)
261 
      real      deltay
 261 
      real      deltay
262 
      parameter (deltay =1.11e5)
 262 
      parameter (deltay =1.11e5)
263 
 
 263 
 
264 
      dvdx = (v1-v0) / ( 2. * deltay * cos(pi180 * lat) )
 264 
      dvdx = (v1-v0) / ( 2. * deltay * cos(pi180 * lat) )
265 
 
 265 
 
266 
      end
 266 
      end
267 
 
 267 
 
268 
c     -------------------------------------------------------------------------
 268 
c     -------------------------------------------------------------------------
269 
c     Zonal derivative of T (DTDX)
 269 
c     Zonal derivative of T (DTDX)
270 
c     -------------------------------------------------------------------------
 270 
c     -------------------------------------------------------------------------
271 
 
 271 
 
272 
      subroutine calc_DTDX (dtdx,t1,t0,lat)
 272 
      subroutine calc_DTDX (dtdx,t1,t0,lat)
273 
 
 273 
 
274 
      implicit none
 274 
      implicit none
275
275
276 
c     Argument declaration
 276 
c     Argument declaration
277 
      real  dtdx   ! Derivative of T in zonal direction [K/m]
 277 
      real  dtdx   ! Derivative of T in zonal direction [K/m]
278 
      real  t1     ! T @ LON + 1 DLON [K]
 278 
      real  t1     ! T @ LON + 1 DLON [K]
279 
      real  t0     ! T @ LON - 1 DLON [K]
 279 
      real  t0     ! T @ LON - 1 DLON [K]
280 
      real  lat    ! Latitude [deg]
 280 
      real  lat    ! Latitude [deg]
281 
 
 281 
 
282 
c     Physical parameters
 282 
c     Physical parameters
283 
      real      pi180
283 
      real      pi180
284 
      parameter (pi180=3.14159/180.)
 284 
      parameter (pi180=3.14159/180.)
285 
      real      deltay
 285 
      real      deltay
286 
      parameter (deltay =1.11e5)
 286 
      parameter (deltay =1.11e5)
287
287
288 
      dtdx = (t1-t0) / ( 2. * deltay * cos(pi180 * lat) )
 288 
      dtdx = (t1-t0) / ( 2. * deltay * cos(pi180 * lat) )
289 
 
 289 
 
290 
      end
 290 
      end
291 
 
 291 
 
292 
c     -------------------------------------------------------------------------
 292 
c     -------------------------------------------------------------------------
293 
c     Zonal derivative of TH (DTHDX)
 293 
c     Zonal derivative of TH (DTHDX)
294 
c     -------------------------------------------------------------------------
 294 
c     -------------------------------------------------------------------------
295 
 
 295 
 
296 
      subroutine calc_DTHDX (dthdx,t1,t0,p,lat)
 296 
      subroutine calc_DTHDX (dthdx,t1,t0,p,lat)
297 
 
 297 
 
298 
      implicit none
 298 
      implicit none
299
299
300 
c     Argument declaration
 300 
c     Argument declaration
301 
      real  dthdx  ! Derivative of TH in zonal direction [K/m]
 301 
      real  dthdx  ! Derivative of TH in zonal direction [K/m]
302 
      real  t1     ! T @ LON + 1 DLON [K]
 302 
      real  t1     ! T @ LON + 1 DLON [K]
303 
      real  t0     ! T @ LON - 1 DLON [K]
 303 
      real  t0     ! T @ LON - 1 DLON [K]
304 
      real  p      ! P [hPa]
 304 
      real  p      ! P [hPa]
305 
      real  lat    ! Latitude [deg]
 305 
      real  lat    ! Latitude [deg]
306 
 
 306 
 
307 
c     Physical parameters
 307 
c     Physical parameters
308 
      real      pi180
308 
      real      pi180
309 
      parameter (pi180=3.14159/180.)
 309 
      parameter (pi180=3.14159/180.)
310 
      real      deltay
 310 
      real      deltay
311 
      parameter (deltay =1.11e5)
 311 
      parameter (deltay =1.11e5)
312 
      real      rdcp,pref
 312 
      real      rdcp,pref
313 
      data      rdcp,pref /0.286,1000./
 313 
      data      rdcp,pref /0.286,1000./
314 
 
 314 
 
315 
      dthdx = (pref/p)**rdcp *
315 
      dthdx = (pref/p)**rdcp *
316 
     >        (t1-t0) / ( 2. * deltay * cos(pi180 * lat) )
 316 
     >        (t1-t0) / ( 2. * deltay * cos(pi180 * lat) )
317 
 
 317 
 
318 
      end
 318 
      end
319 
 
 319 
 
320 
c     -------------------------------------------------------------------------
 320 
c     -------------------------------------------------------------------------
321 
c     Meridional derivative of U (DUDY)
 321 
c     Meridional derivative of U (DUDY)
322 
c     -------------------------------------------------------------------------
 322 
c     -------------------------------------------------------------------------
323 
 
 323 
 
324 
      subroutine calc_DUDY (dudy,u1,u0)
 324 
      subroutine calc_DUDY (dudy,u1,u0)
325 
 
 325 
 
326 
      implicit none
 326 
      implicit none
327
327
328 
c     Argument declaration
 328 
c     Argument declaration
329 
      real  dudy   ! Derivative of U in meridional direction [s^-1]
 329 
      real  dudy   ! Derivative of U in meridional direction [s^-1]
330 
      real  u1     ! U @ LAT + 1 DLAT [m/s]
 330 
      real  u1     ! U @ LAT + 1 DLAT [m/s]
331 
      real  u0     ! U @ LAT - 1 DLAT [m/s]
 331 
      real  u0     ! U @ LAT - 1 DLAT [m/s]
332 
 
 332 
 
333 
c     Physical parameters
 333 
c     Physical parameters
334 
      real      deltay
 334 
      real      deltay
335 
      parameter (deltay =1.11e5)
 335 
      parameter (deltay =1.11e5)
336
336
337 
      dudy = (u1-u0) / ( 2. * deltay )
 337 
      dudy = (u1-u0) / ( 2. * deltay )
338 
 
 338 
 
339 
      end
 339 
      end
340 
 
 340 
 
341 
c     -------------------------------------------------------------------------
 341 
c     -------------------------------------------------------------------------
342 
c     Meridional derivative of V (DVDY)
 342 
c     Meridional derivative of V (DVDY)
343 
c     -------------------------------------------------------------------------
 343 
c     -------------------------------------------------------------------------
344 
 
 344 
 
345 
      subroutine calc_DVDY (dvdy,v1,v0)
 345 
      subroutine calc_DVDY (dvdy,v1,v0)
346 
 
 346 
 
347 
      implicit none
 347 
      implicit none
348
348
349 
c     Argument declaration
 349 
c     Argument declaration
350 
      real  dvdy   ! Derivative of V in meridional direction [s^-1]
 350 
      real  dvdy   ! Derivative of V in meridional direction [s^-1]
351 
      real  v1     ! V @ LAT + 1 DLAT [m/s]
 351 
      real  v1     ! V @ LAT + 1 DLAT [m/s]
352 
      real  v0     ! V @ LAT - 1 DLAT [m/s]
 352 
      real  v0     ! V @ LAT - 1 DLAT [m/s]
353 
 
 353 
 
354 
c     Physical parameters
 354 
c     Physical parameters
355 
      real      deltay
 355 
      real      deltay
356 
      parameter (deltay =1.11e5)
 356 
      parameter (deltay =1.11e5)
357
357
358 
      dvdy = (v1-v0) / ( 2. * deltay )
 358 
      dvdy = (v1-v0) / ( 2. * deltay )
359 
 
 359 
 
360 
      end
 360 
      end
361 
 
 361 
 
362 
c     -------------------------------------------------------------------------
 362 
c     -------------------------------------------------------------------------
363 
c     Meridional derivative of T (DTDY)
 363 
c     Meridional derivative of T (DTDY)
364 
c     -------------------------------------------------------------------------
 364 
c     -------------------------------------------------------------------------
365 
 
 365 
 
366 
      subroutine calc_DTDY (dtdy,t1,t0)
 366 
      subroutine calc_DTDY (dtdy,t1,t0)
367 
 
 367 
 
368 
      implicit none
 368 
      implicit none
369
369
370 
c     Argument declaration
 370 
c     Argument declaration
371 
      real  dtdy   ! Derivative of T in meridional direction [K/m]
 371 
      real  dtdy   ! Derivative of T in meridional direction [K/m]
372 
      real  t1     ! T @ LAT + 1 DLAT [K]
 372 
      real  t1     ! T @ LAT + 1 DLAT [K]
373 
      real  t0     ! T @ LAT - 1 DLAT [K]
 373 
      real  t0     ! T @ LAT - 1 DLAT [K]
374 
 
 374 
 
375 
c     Physical parameters
 375 
c     Physical parameters
376 
      real      deltay
 376 
      real      deltay
377 
      parameter (deltay =1.11e5)
 377 
      parameter (deltay =1.11e5)
378
378
379 
      dtdy = (t1-t0) / ( 2. * deltay )
 379 
      dtdy = (t1-t0) / ( 2. * deltay )
380 
 
 380 
 
381 
      end
 381 
      end
382 
 
 382 
 
383 
c     -------------------------------------------------------------------------
 383 
c     -------------------------------------------------------------------------
384 
c     Meridional derivative of TH (DTHDY)
 384 
c     Meridional derivative of TH (DTHDY)
385 
c     -------------------------------------------------------------------------
 385 
c     -------------------------------------------------------------------------
386 
 
 386 
 
387 
      subroutine calc_DTHDY (dthdy,t1,t0,p)
 387 
      subroutine calc_DTHDY (dthdy,t1,t0,p)
388 
 
 388 
 
389 
      implicit none
 389 
      implicit none
390
390
391 
c     Argument declaration
 391 
c     Argument declaration
392 
      real  dthdy  ! Derivative of TH in meridional direction [K/m]
 392 
      real  dthdy  ! Derivative of TH in meridional direction [K/m]
393 
      real  t1     ! TH @ LAT + 1 DLAT [K]
 393 
      real  t1     ! TH @ LAT + 1 DLAT [K]
394 
      real  t0     ! TH @ LAT - 1 DLAT [K]
 394 
      real  t0     ! TH @ LAT - 1 DLAT [K]
395 
      real  p      ! P [hPa]
 395 
      real  p      ! P [hPa]
396 
 
 396 
 
397 
c     Physical parameters
 397 
c     Physical parameters
398 
      real      deltay
 398 
      real      deltay
399 
      parameter (deltay =1.11e5)
 399 
      parameter (deltay =1.11e5)
400 
      real      rdcp,pref
 400 
      real      rdcp,pref
401 
      data      rdcp,pref /0.286,1000./
 401 
      data      rdcp,pref /0.286,1000./
402 
 
 402 
 
403 
      dthdy = (pref/p)**rdcp * (t1-t0) / ( 2. * deltay )
 403 
      dthdy = (pref/p)**rdcp * (t1-t0) / ( 2. * deltay )
404 
 
 404 
 
405 
      end
 405 
      end
406 
 
 406 
 
407 
c     -------------------------------------------------------------------------
 407 
c     -------------------------------------------------------------------------
408 
c     Wind shear of U (DUDP)
 408 
c     Wind shear of U (DUDP)
409 
c     -------------------------------------------------------------------------
 409 
c     -------------------------------------------------------------------------
410 
 
 410 
 
411 
      subroutine calc_DUDP (dudp,u1,u0,p1,p0)
 411 
      subroutine calc_DUDP (dudp,u1,u0,p1,p0)
412 
 
 412 
 
413 
      implicit none
 413 
      implicit none
414
414
415 
c     Argument declaration
 415 
c     Argument declaration
416 
      real  dudp   ! Wind shear [m/s per Pa]
 416 
      real  dudp   ! Wind shear [m/s per Pa]
417 
      real  u1     ! U @ P + 1 DP [m/s]
 417 
      real  u1     ! U @ P + 1 DP [m/s]
418 
      real  u0     ! U @ P - 1 DP [m/s]
 418 
      real  u0     ! U @ P - 1 DP [m/s]
419 
      real  p1     ! P + 1 DP [hPa]
 419 
      real  p1     ! P + 1 DP [hPa]
420 
      real  p0     ! P - 1 DP [hPa]
 420 
      real  p0     ! P - 1 DP [hPa]
421 
 
 421 
 
422 
      dudp = 0.01 * (u1-u0) / (p1-p0)
 422 
      dudp = 0.01 * (u1-u0) / (p1-p0)
423 
 
 423 
 
424 
      end
 424 
      end
425 
 
 425 
 
426 
c     -------------------------------------------------------------------------
 426 
c     -------------------------------------------------------------------------
427 
c     Wind shear of V (DVDP)
 427 
c     Wind shear of V (DVDP)
428 
c     -------------------------------------------------------------------------
 428 
c     -------------------------------------------------------------------------
429 
 
 429 
 
430 
      subroutine calc_DVDP (dvdp,v1,v0,p1,p0)
 430 
      subroutine calc_DVDP (dvdp,v1,v0,p1,p0)
431 
 
 431 
 
432 
      implicit none
 432 
      implicit none
433 
 
 433 
 
434 
c     Argument declaration
 434 
c     Argument declaration
435 
      real  dvdp   ! Wind shear [m/s per Pa]
 435 
      real  dvdp   ! Wind shear [m/s per Pa]
436 
      real  v1     ! V @ P + 1 DP [m/s]
 436 
      real  v1     ! V @ P + 1 DP [m/s]
437 
      real  v0     ! V @ P - 1 DP [m/s]
 437 
      real  v0     ! V @ P - 1 DP [m/s]
438 
      real  p1     ! P + 1 DP [hPa]
 438 
      real  p1     ! P + 1 DP [hPa]
439 
      real  p0     ! P - 1 DP [hPa]
 439 
      real  p0     ! P - 1 DP [hPa]
440 
 
 440 
 
441 
      dvdp = 0.01 * (v1-v0) / (p1-p0)
 441 
      dvdp = 0.01 * (v1-v0) / (p1-p0)
442 
 
 442 
 
443 
      end
 443 
      end
444 
 
 444 
 
445 
c     -------------------------------------------------------------------------
 445 
c     -------------------------------------------------------------------------
446 
c     Vertical derivative of T (DTDP)
 446 
c     Vertical derivative of T (DTDP)
447 
c     -------------------------------------------------------------------------
 447 
c     -------------------------------------------------------------------------
448 
 
 448 
 
449 
      subroutine calc_DTDP (dtdp,t1,t0,p1,p0)
 449 
      subroutine calc_DTDP (dtdp,t1,t0,p1,p0)
450
450
451 
      implicit none
 451 
      implicit none
452 
 
 452 
 
453 
c     Argument declaration
 453 
c     Argument declaration
454 
      real  dtdp   ! Vertical derivative of T [K/Pa]
 454 
      real  dtdp   ! Vertical derivative of T [K/Pa]
455 
      real  t1     ! T @ P + 1 DP [K]
 455 
      real  t1     ! T @ P + 1 DP [K]
456 
      real  t0     ! T @ P - 1 DP [K]
 456 
      real  t0     ! T @ P - 1 DP [K]
457 
      real  p1     ! P + 1 DP [hPa]
 457 
      real  p1     ! P + 1 DP [hPa]
458 
      real  p0     ! P - 1 DP [hPa]
 458 
      real  p0     ! P - 1 DP [hPa]
459 
 
 459 
 
460 
      dtdp = 0.01 * (t1-t0) / (p1-p0)
 460 
      dtdp = 0.01 * (t1-t0) / (p1-p0)
461 
 
 461 
 
462 
      end
 462 
      end
463 
 
 463 
 
464 
c     -------------------------------------------------------------------------
 464 
c     -------------------------------------------------------------------------
465 
c     Vertical derivative of TH (DTHDP)
 465 
c     Vertical derivative of TH (DTHDP)
466 
c     -------------------------------------------------------------------------
 466 
c     -------------------------------------------------------------------------
467 
 
 467 
 
468 
      subroutine calc_DTHDP (dthdp,t1,t0,p1,p0,p,t)
 468 
      subroutine calc_DTHDP (dthdp,t1,t0,p1,p0,p,t)
469 
 
 469 
 
470 
      implicit none
 470 
      implicit none
471
471
472 
c     Argument declaration
 472 
c     Argument declaration
473 
      real  dthdp  ! Vertical derivative of TH [K/Pa]
 473 
      real  dthdp  ! Vertical derivative of TH [K/Pa]
474 
      real  t1     ! T @ P + 1 DP [K]
 474 
      real  t1     ! T @ P + 1 DP [K]
475 
      real  t0     ! T @ P - 1 DP [K]
 475 
      real  t0     ! T @ P - 1 DP [K]
476 
      real  t      ! T [K]
 476 
      real  t      ! T [K]
477 
      real  p1     ! P + 1 DP [hPa]
 477 
      real  p1     ! P + 1 DP [hPa]
478 
      real  p0     ! P - 1 DP [hPa]
 478 
      real  p0     ! P - 1 DP [hPa]
479 
      real  p      ! P [hPa]
 479 
      real  p      ! P [hPa]
480 
 
 480 
 
481 
c     Physical parameters
 481 
c     Physical parameters
482 
      real      rdcp,tzero,pref
 482 
      real      rdcp,tzero,pref
483 
      data      rdcp,tzero,pref /0.286,273.15,1000./
 483 
      data      rdcp,tzero,pref /0.286,273.15,1000./
484
484
485 
c     Auxiliary variables
485 
c     Auxiliary variables
486 
      real      tk1,tk0,tk
 486 
      real      tk1,tk0,tk
487 
 
 487 
 
488 
      if (t0.lt.100.) then
 488 
      if (t0.lt.100.) then
489 
         tk0 = t0 + tzero
 489 
         tk0 = t0 + tzero
490 
      endif
 490 
      endif
491 
      if (t1.lt.100.) then
 491 
      if (t1.lt.100.) then
492 
         tk1 = t1 + tzero
 492 
         tk1 = t1 + tzero
493 
      endif
493 
      endif
494 
      if (t.lt.100.) then
 494 
      if (t.lt.100.) then
495 
         tk  = t  + tzero
 495 
         tk  = t  + tzero
496 
      endif
496 
      endif
497 
 
 497 
 
498 
      dthdp = 0.01*(pref/p)**rdcp *
498 
      dthdp = 0.01*(pref/p)**rdcp *
499 
     >      ( (tk1-tk0)/(p1-p0) - rdcp * tk/p )
499 
     >      ( (tk1-tk0)/(p1-p0) - rdcp * tk/p )
500 
 
 500 
 
501 
      end
 501 
      end
502 
 
 502 
 
503 
c     -------------------------------------------------------------------------
 503 
c     -------------------------------------------------------------------------
504 
c     Squared Brunt-Vaisäla frequency (NSQ)
 504 
c     Squared Brunt-Vaisäla frequency (NSQ)
505 
c     -------------------------------------------------------------------------
 505 
c     -------------------------------------------------------------------------
506 
 
 506 
 
507 
      subroutine calc_NSQ (nsq,dthdp,th,rho)
 507 
      subroutine calc_NSQ (nsq,dthdp,th,rho)
508 
 
 508 
 
509 
      implicit none
 509 
      implicit none
510
510
511 
c     Argument declaration
 511 
c     Argument declaration
512 
      real  nsq    ! Squared Brunt-Vaisäla frequency [s^-1]
 512 
      real  nsq    ! Squared Brunt-Vaisäla frequency [s^-1]
513 
      real  dthdp  ! D(TH)/DP [K/Pa]
 513 
      real  dthdp  ! D(TH)/DP [K/Pa]
514 
      real  th     ! K
 514 
      real  th     ! K
515 
      real  rho    ! Density [kg m^-3]
515 
      real  rho    ! Density [kg m^-3]
516 
 
 516 
 
517 
c     Physical parameters
 517 
c     Physical parameters
518 
      real      g
518 
      real      g
519 
      parameter (g=9.81)
 519 
      parameter (g=9.81)
520 
 
 520 
 
521 
      nsq = -g**2/th * rho * dthdp
 521 
      nsq = -g**2/th * rho * dthdp
522 
 
 522 
 
523 
      end
 523 
      end
524 
 
 524 
 
525 
c     -------------------------------------------------------------------------
 525 
c     -------------------------------------------------------------------------
526 
c     Relative vorticity (RELVORT)
 526 
c     Relative vorticity (RELVORT)
527 
c     -------------------------------------------------------------------------
 527 
c     -------------------------------------------------------------------------
528 
 
 528 
 
529 
      subroutine calc_RELVORT (relvort,dudy,dvdx,u,lat)
 529 
      subroutine calc_RELVORT (relvort,dudy,dvdx,u,lat)
530 
 
 530 
 
531 
      implicit none
 531 
      implicit none
532
532
533 
c     Argument declaration
 533 
c     Argument declaration
534 
      real  relvort    ! Relative vorticity [s^-1]
 534 
      real  relvort    ! Relative vorticity [s^-1]
535 
      real  u          ! Zonal wind component [m/s]
 535 
      real  u          ! Zonal wind component [m/s]
536 
      real  dudy       ! du/dy [s^-1]
 536 
      real  dudy       ! du/dy [s^-1]
537 
      real  dvdx       ! dv/dx [s^-1]
 537 
      real  dvdx       ! dv/dx [s^-1]
538 
      real  lat        ! Latitude [deg]
 538 
      real  lat        ! Latitude [deg]
539 
 
 539 
 
540 
c     Physical parameters
 540 
c     Physical parameters
541 
      real      pi180
 541 
      real      pi180
542 
      parameter (pi180=3.14159/180.)
 542 
      parameter (pi180=3.14159/180.)
543 
      real      deltay
 543 
      real      deltay
544 
      parameter (deltay =1.11e5)
 544 
      parameter (deltay =1.11e5)
545 
 
 545 
 
546 
      relvort = dvdx - dudy + u * pi180/deltay * tan(pi180 * lat)
546 
      relvort = dvdx - dudy + u * pi180/deltay * tan(pi180 * lat)
547
547
548 
      end
 548 
      end
549 
 
 549 
 
550 
c     -------------------------------------------------------------------------
 550 
c     -------------------------------------------------------------------------
551 
c     Absolute vorticity (ABSVORT)
 551 
c     Absolute vorticity (ABSVORT)
552 
c     -------------------------------------------------------------------------
 552 
c     -------------------------------------------------------------------------
553 
 
 553 
 
554 
      subroutine calc_ABSVORT (absvort,dudy,dvdx,u,lat)
 554 
      subroutine calc_ABSVORT (absvort,dudy,dvdx,u,lat)
555 
 
 555 
 
556 
      implicit none
 556 
      implicit none
557 
 
 557 
 
558 
c     Argument declaration
 558 
c     Argument declaration
559 
      real  absvort    ! Absolute vorticity [s^-1]
 559 
      real  absvort    ! Absolute vorticity [s^-1]
560 
      real  u          ! Zonal wind component [m/s]
 560 
      real  u          ! Zonal wind component [m/s]
561 
      real  dudy       ! du/dy [s^-1]
 561 
      real  dudy       ! du/dy [s^-1]
562 
      real  dvdx       ! dv/dx [s^-1]
 562 
      real  dvdx       ! dv/dx [s^-1]
563 
      real  lat        ! Latitude [deg]
 563 
      real  lat        ! Latitude [deg]
564 
 
 564 
 
565 
c     Physical parameters
 565 
c     Physical parameters
566 
      real      pi180
 566 
      real      pi180
567 
      parameter (pi180=3.14159/180.)
 567 
      parameter (pi180=3.14159/180.)
568 
      real      deltay
 568 
      real      deltay
569 
      parameter (deltay =1.11e5)
 569 
      parameter (deltay =1.11e5)
570 
      real      omega
 570 
      real      omega
571 
      parameter (omega=7.292e-5)
 571 
      parameter (omega=7.292e-5)
572 
 
 572 
 
573 
      absvort = dvdx - dudy + u * pi180/deltay * tan(pi180 * lat) +
 573 
      absvort = dvdx - dudy + u * pi180/deltay * tan(pi180 * lat) +
574 
     >          2. * omega * sin(pi180 * lat)
574 
     >          2. * omega * sin(pi180 * lat)
575
575
576 
      end
 576 
      end
577 
 
 577 
 
578 
c     -------------------------------------------------------------------------
 578 
c     -------------------------------------------------------------------------
579 
c     Divergence (DIV)
 579 
c     Divergence (DIV)
580 
c     -------------------------------------------------------------------------
 580 
c     -------------------------------------------------------------------------
581 
 
 581 
 
582 
      subroutine calc_DIV (div,dudx,dvdy,v,lat)
 582 
      subroutine calc_DIV (div,dudx,dvdy,v,lat)
583 
 
 583 
 
584 
      implicit none
 584 
      implicit none
585
585
586 
c     Argument declaration
 586 
c     Argument declaration
587 
      real  div        ! Divergence [s^-1]
 587 
      real  div        ! Divergence [s^-1]
588 
      real  v          ! Meridional wind component [m/s]
 588 
      real  v          ! Meridional wind component [m/s]
589 
      real  dudx       ! du/dx [s^-1]
 589 
      real  dudx       ! du/dx [s^-1]
590 
      real  dvdy       ! dv/dy [s^-1]
 590 
      real  dvdy       ! dv/dy [s^-1]
591 
      real  lat        ! Latitude [deg]
 591 
      real  lat        ! Latitude [deg]
592 
 
 592 
 
593 
c     Physical parameters
 593 
c     Physical parameters
594 
      real      pi180
 594 
      real      pi180
595 
      parameter (pi180=3.14159/180.)
 595 
      parameter (pi180=3.14159/180.)
596 
      real      deltay
 596 
      real      deltay
597 
      parameter (deltay =1.11e5)
 597 
      parameter (deltay =1.11e5)
598 
 
 598 
 
599 
      div = dudx + dvdy - v * pi180/deltay * tan(pi180 * lat)
599 
      div = dudx + dvdy - v * pi180/deltay * tan(pi180 * lat)
600
600
601 
      end
 601 
      end
602 
 
 602 
 
603 
c     -------------------------------------------------------------------------
 603 
c     -------------------------------------------------------------------------
604 
c     Deformation (DEF)
 604 
c     Deformation (DEF)
605 
c     -------------------------------------------------------------------------
 605 
c     -------------------------------------------------------------------------
606 
 
 606 
 
607 
      subroutine calc_DEF (def,dudx,dvdx,dudy,dvdy)
 607 
      subroutine calc_DEF (def,dudx,dvdx,dudy,dvdy)
608 
 
 608 
 
609 
      implicit none
 609 
      implicit none
610 
 
 610 
 
611 
c     Argument declaration
 611 
c     Argument declaration
612 
      real  def        ! Deformation [s^-1]
 612 
      real  def        ! Deformation [s^-1]
613 
      real  dudx       ! du/dx [s^-1]
 613 
      real  dudx       ! du/dx [s^-1]
614 
      real  dvdx       ! dv/dy [s^-1]
 614 
      real  dvdx       ! dv/dy [s^-1]
615 
      real  dudy       ! du/dx [s^-1]
 615 
      real  dudy       ! du/dx [s^-1]
616 
      real  dvdy       ! dv/dy [s^-1]
 616 
      real  dvdy       ! dv/dy [s^-1]
617 
 
 617 
 
618 
c     Physical parameters
 618 
c     Physical parameters
619 
      real      pi180
 619 
      real      pi180
620 
      parameter (pi180=3.14159/180.)
 620 
      parameter (pi180=3.14159/180.)
621 
 
 621 
 
622 
      def = sqrt( (dvdx+dudy)**2 + (dudx-dvdy)**2 )
 622 
      def = sqrt( (dvdx+dudy)**2 + (dudx-dvdy)**2 )
623
623
624 
      end
 624 
      end
625 
 
 625 
 
626 
c     -------------------------------------------------------------------------
 626 
c     -------------------------------------------------------------------------
627 
c     Potential Vorticity (PV)
 627 
c     Potential Vorticity (PV)
628 
c     -------------------------------------------------------------------------
 628 
c     -------------------------------------------------------------------------
629 
 
 629 
 
630 
      subroutine calc_PV (pv,absvort,dthdp,dudp,dvdp,dthdx,dthdy)
 630 
      subroutine calc_PV (pv,absvort,dthdp,dudp,dvdp,dthdx,dthdy)
631 
 
 631 
 
632 
      implicit none
 632 
      implicit none
633 
 
 633 
 
634 
c     Argument declaration
 634 
c     Argument declaration
635 
      real  pv         ! Ertel-PV [PVU]
 635 
      real  pv         ! Ertel-PV [PVU]
636 
      real  absvort    ! Absolute vorticity [s^-1]
 636 
      real  absvort    ! Absolute vorticity [s^-1]
637 
      real  dthdp      ! dth/dp [K/Pa]
 637 
      real  dthdp      ! dth/dp [K/Pa]
638 
      real  dudp       ! du/dp [m/s per Pa]
 638 
      real  dudp       ! du/dp [m/s per Pa]
639 
      real  dvdp       ! dv/dp [m/s per Pa]
 639 
      real  dvdp       ! dv/dp [m/s per Pa]
640 
      real  dthdx      ! dth/dx [K/m]
 640 
      real  dthdx      ! dth/dx [K/m]
641 
      real  dthdy      ! dth/dy [K/m]
 641 
      real  dthdy      ! dth/dy [K/m]
642 
 
 642 
 
643 
c     Physical and numerical parameters
 643 
c     Physical and numerical parameters
644 
      real      scale
 644 
      real      scale
645 
      parameter (scale=1.E6)
 645 
      parameter (scale=1.E6)
646 
      real      g
 646 
      real      g
647 
      parameter (g=9.80665)
 647 
      parameter (g=9.80665)
648 
 
 648 
 
649 
      pv = -scale * g * ( absvort * dthdp + dudp * dthdy - dvdp * dthdx)
 649 
      pv = -scale * g * ( absvort * dthdp + dudp * dthdy - dvdp * dthdx)
650 
 
 650 
 
651 
      end
 651 
      end
652 
 
 652 
 
653 
c     -------------------------------------------------------------------------
 653 
c     -------------------------------------------------------------------------
654 
c     Richardson number (RI)
 654 
c     Richardson number (RI)
655 
c     -------------------------------------------------------------------------
 655 
c     -------------------------------------------------------------------------
656 
 
 656 
 
657 
      subroutine calc_RI (ri,dudp,dvdp,nsq,rho)
 657 
      subroutine calc_RI (ri,dudp,dvdp,nsq,rho)
658 
 
 658 
 
659 
      implicit none
 659 
      implicit none
660 
 
 660 
 
661 
c     Argument declaration
 661 
c     Argument declaration
662 
      real  ri         ! Richardson number
 662 
      real  ri         ! Richardson number
663 
      real  dudp       ! Du/Dp [m/s per Pa]
 663 
      real  dudp       ! Du/Dp [m/s per Pa]
664 
      real  dvdp       ! Dv/Dp [m/s per Pa]
 664 
      real  dvdp       ! Dv/Dp [m/s per Pa]
665 
      real  nsq        ! Squared Brunt-Vailälä frequency [s^-1]
 665 
      real  nsq        ! Squared Brunt-Vailälä frequency [s^-1]
666 
      real  rho        ! Density [kg/m^3]
 666 
      real  rho        ! Density [kg/m^3]
667 
 
 667 
 
668 
c     Physical and numerical parameters
 668 
c     Physical and numerical parameters
669 
      real      g
 669 
      real      g
670 
      parameter (g=9.80665)
 670 
      parameter (g=9.80665)
671 
 
 671 
 
672 
      ri = nsq / ( dudp**2 + dvdp**2 ) / ( rho * g )**2
 672 
      ri = nsq / ( dudp**2 + dvdp**2 ) / ( rho * g )**2
673 
 
 673 
 
674 
      end
 674 
      end
675 
 
 675 
 
676 
c     -------------------------------------------------------------------------
 676 
c     -------------------------------------------------------------------------
677 
c     Ellrod and Knapp's turbulence indicator (TI)
 677 
c     Ellrod and Knapp's turbulence indicator (TI)
678 
c     -------------------------------------------------------------------------
 678 
c     -------------------------------------------------------------------------
679 
 
 679 
 
680 
      subroutine calc_TI (ti,def,dudp,dvdp,rho)
 680 
      subroutine calc_TI (ti,def,dudp,dvdp,rho)
681
681
682 
      implicit none
 682 
      implicit none
683 
 
 683 
 
684 
c     Argument declaration
 684 
c     Argument declaration
685 
      real  ti         ! Turbulence idicator
 685 
      real  ti         ! Turbulence idicator
686 
      real  def        ! Deformation [s^-1]
 686 
      real  def        ! Deformation [s^-1]
687 
      real  dudp       ! Du/Dp [m/s per Pa]
 687 
      real  dudp       ! Du/Dp [m/s per Pa]
688 
      real  dvdp       ! Dv/Dp [m/s per Pa]
 688 
      real  dvdp       ! Dv/Dp [m/s per Pa]
689 
      real  rho        ! Density [kg/m^3]
 689 
      real  rho        ! Density [kg/m^3]
690 
 
 690 
 
691 
c     Physical and numerical parameters
 691 
c     Physical and numerical parameters
692 
      real      g
 692 
      real      g
693 
      parameter (g=9.80665)
 693 
      parameter (g=9.80665)
694 
 
 694 
 
695 
      ti = def * sqrt ( dudp**2 + dvdp**2 ) * ( rho * g )
 695 
      ti = def * sqrt ( dudp**2 + dvdp**2 ) * ( rho * g )
696 
 
 696 
 
697 
      end
 697 
      end
698 
 
 698 
 
699 
c     -------------------------------------------------------------------------
 699 
c     -------------------------------------------------------------------------
700 
c     Distance from starting position
 700 
c     Distance from starting position
701 
c     -------------------------------------------------------------------------
 701 
c     -------------------------------------------------------------------------
702 
 
 702 
 
703 
      subroutine calc_DIST0 (dist0,lon0,lat0,lon1,lat1)
 703 
      subroutine calc_DIST0 (dist0,lon0,lat0,lon1,lat1)
704 
 
 704 
 
705 
      implicit none
 705 
      implicit none
706 
 
 706 
 
707 
c     Argument declaration
 707 
c     Argument declaration
708 
      real  dist0      ! Distance from starting position [km]
 708 
      real  dist0      ! Distance from starting position [km]
709 
      real  lon0,lat0  ! Starting position
 709 
      real  lon0,lat0  ! Starting position
710 
      real  lon1,lat1  ! New position
710 
      real  lon1,lat1  ! New position
711 
 
 711 
 
712 
c     Externals
712 
c     Externals
713 
      real     sdis
 713 
      real     sdis
714 
      external sdis
 714 
      external sdis
715 
 
 715 
 
716 
      dist0 = sdis(lon0,lat0,lon1,lat1)
 716 
      dist0 = sdis(lon0,lat0,lon1,lat1)
717
717
718 
      end
 718 
      end
719 
 
 719 
 
720 
c     -------------------------------------------------------------------------
 720 
c     -------------------------------------------------------------------------
721 
c     Heading of the trajectory (HEAD)
 721 
c     Heading of the trajectory (HEAD)
722 
c     -------------------------------------------------------------------------
 722 
c     -------------------------------------------------------------------------
723 
 
 723 
 
724 
      subroutine calc_HEAD (head,lon0,lat0,lon1,lat1)
 724 
      subroutine calc_HEAD (head,lon0,lat0,lon1,lat1)
725
725
726 
      implicit none
 726 
      implicit none
727 
 
 727 
 
728 
c     Argument declaration
 728 
c     Argument declaration
729 
      real  head       ! Heading angle (in deg) relativ to zonal direction
 729 
      real  head       ! Heading angle (in deg) relativ to zonal direction
730 
      real  lon0,lat0  ! Starting position
 730 
      real  lon0,lat0  ! Starting position
731 
      real  lon1,lat1  ! New position
731 
      real  lon1,lat1  ! New position
732 
 
 732 
 
733 
c     Physical parameters
 733 
c     Physical parameters
734 
      real      pi180
 734 
      real      pi180
735 
      parameter (pi180=3.14159/180.)
 735 
      parameter (pi180=3.14159/180.)
736 
 
 736 
 
737 
c     Auixiliary variables
 737 
c     Auixiliary variables
738 
      real     dx,dy
 738 
      real     dx,dy
- 
 
 739 
      real     dlon
739 
 
 740 
 
- 
 
 741 
      dlon = lon1-lon0
- 
 
 742 
      if ( dlon.gt.180.  ) dlon = dlon - 360.
- 
 
 743 
      if ( dlon.lt.-180. ) dlon = dlon + 360.
- 
 
 744
740 
      dx = (lon1-lon0) * cos(pi180*0.5*(lat0+lat1))
 745 
      dx = dlon * cos(pi180*0.5*(lat0+lat1))
741 
      dy = lat1-lat0
 746 
      dy = lat1-lat0
742 
 
 747 
 
743 
      call getangle(1.,0.,dx,dy,head)
 748 
      call getangle(1.,0.,dx,dy,head)
744 
 
 749 
 
- 
 
 750 
      end
- 
 
 751
- 
 
 752 
c     -------------------------------------------------------------------------
- 
 
 753 
c     Directional change of the trajectory (HEAD)
- 
 
 754 
c     -------------------------------------------------------------------------
- 
 
 755 
 
- 
 
 756 
      subroutine calc_DANGLE (dangle,lon0,lat0,lon1,lat1,lon2,lat2)
- 
 
 757
- 
 
 758 
      implicit none
- 
 
 759 
 
- 
 
 760 
c     Argument declaration
- 
 
 761 
      real  dangle     ! Directiopanl change
- 
 
 762 
      real  lon0,lat0  ! t-1
- 
 
 763 
      real  lon1,lat1  ! t
- 
 
 764 
      real  lon2,lat2  ! t+1
- 
 
 765
- 
 
 766 
c     Physical parameters
- 
 
 767 
      real      pi180
- 
 
 768 
      parameter (pi180=3.14159/180.)
- 
 
 769 
      real      eps
- 
 
 770 
      parameter (eps=0.00001)
- 
 
 771 
 
- 
 
 772 
c     Auixiliary variables
- 
 
 773 
      real     dx1,dy1,dx2,dy2,norm,cross,dlon1,dlon2
- 
 
 774
- 
 
 775 
      dlon1 = lon1 - lon0
- 
 
 776 
      if ( dlon1.gt.180.  ) dlon1 = dlon1 - 360.
- 
 
 777 
      if ( dlon1.lt.-180. ) dlon1 = dlon1 + 360.
- 
 
 778 
      dlon2 = lon2 - lon1
- 
 
 779 
      if ( dlon2.gt.180.  ) dlon2 = dlon2 - 360.
- 
 
 780 
      if ( dlon2.lt.-180. ) dlon2 = dlon2 + 360.
- 
 
 781
- 
 
 782 
      dx1 = dlon1 * cos(pi180*0.5*(lat1+lat0))
- 
 
 783 
      dy1 = lat1 - lat0
- 
 
 784 
      dx2 = dlon2 * cos(pi180*0.5*(lat2+lat1))
- 
 
 785 
      dy2 = lat2 - lat1
- 
 
 786
- 
 
 787 
      norm = sqrt( (dx1**2 + dy1**2) * (dx2**2 + dy2**2) )
- 
 
 788 
 
- 
 
 789 
      if ( norm.gt.eps ) then
- 
 
 790 
         cross = ( dx1 * dy2 - dy1 * dx2 ) / norm
- 
 
 791 
         if ( cross.ge.1. ) then
- 
 
 792 
            dangle = 90.
- 
 
 793 
         elseif (cross.le.-1.) then
- 
 
 794 
            dangle = -90.
- 
 
 795 
         else
- 
 
 796 
            dangle = 1./pi180 * asin( cross )
- 
 
 797 
         endif
- 
 
 798 
      else
- 
 
 799 
         dangle = -999.
- 
 
 800 
      endif
- 
 
 801 
 
745 
      end
 802 
      end
746 
 
 803 
 
747 
c
804 
c
748 
c     *************************************************************************
 805 
c     *************************************************************************
749 
c     Auxiliary subroutines and functions
 806 
c     Auxiliary subroutines and functions
750 
c     *************************************************************************
 807 
c     *************************************************************************
751
808
752 
c     -------------------------------------------------------------------------
 809 
c     -------------------------------------------------------------------------
753 
c     Saturation vapor pressure over water
 810 
c     Saturation vapor pressure over water
754 
c     -------------------------------------------------------------------------
 811 
c     -------------------------------------------------------------------------
755 
 
 812 
 
756 
      real function esat(t)
 813 
      real function esat(t)
757
814
758 
C     This function returns the saturation vapor pressure over water
 815 
C     This function returns the saturation vapor pressure over water
759 
c     (mb) given the temperature (Kelvin).
 816 
c     (mb) given the temperature (Kelvin).
760 
C     The algorithm is due to Nordquist, W. S. ,1973: "Numerical
 817 
C     The algorithm is due to Nordquist, W. S. ,1973: "Numerical
761 
C     Approximations of Selected Meteorological Parameters for Cloud
 818 
C     Approximations of Selected Meteorological Parameters for Cloud
762 
C     Physics Problems" ECOM-5475, Atmospheric Sciences Laboratory,
 819 
C     Physics Problems" ECOM-5475, Atmospheric Sciences Laboratory,
763 
c     U. S. Army Electronics Command, White Sands Missile Range,
 820 
c     U. S. Army Electronics Command, White Sands Missile Range,
764 
c     New Mexico 88002.
 821 
c     New Mexico 88002.
765
822
766 
      real p1,p2,c1,t
 823 
      real p1,p2,c1,t
767
824
768 
      p1=11.344-0.0303998*t
 825 
      p1=11.344-0.0303998*t
769 
      p2=3.49149-1302.8844/t
 826 
      p2=3.49149-1302.8844/t
770 
      c1=23.832241-5.02808*log10(t)
 827 
      c1=23.832241-5.02808*log10(t)
771 
      esat=10.**(c1-1.3816e-7*10.**p1+8.1328e-3*10.**p2-2949.076/t)
 828 
      esat=10.**(c1-1.3816e-7*10.**p1+8.1328e-3*10.**p2-2949.076/t)
772 
 
 829 
 
773 
      end
 830 
      end
774 
 
 831 
 
775 
c     --------------------------------------------------------------------------
 832 
c     --------------------------------------------------------------------------
776 
c     Angle between two vectors
 833 
c     Angle between two vectors
777 
c     --------------------------------------------------------------------------
 834 
c     --------------------------------------------------------------------------
778 
 
 835 
 
779 
      SUBROUTINE getangle (ux1,uy1,ux2,uy2,angle)
 836 
      SUBROUTINE getangle (ux1,uy1,ux2,uy2,angle)
780 
 
 837 
 
781 
c     Given two vectors <ux1,uy1> and <ux2,uy2>, determine the angle (in deg)
 838 
c     Given two vectors <ux1,uy1> and <ux2,uy2>, determine the angle (in deg)
782 
c     between the two vectors.
 839 
c     between the two vectors.
783 
 
 840 
 
784 
      implicit none
 841 
      implicit none
785 
 
 842 
 
786 
c     Declaration of subroutine parameters
 843 
c     Declaration of subroutine parameters
787 
      real ux1,uy1
 844 
      real ux1,uy1
788 
      real ux2,uy2
 845 
      real ux2,uy2
789 
      real angle
 846 
      real angle
790 
 
 847 
 
791 
c     Auxiliary variables and parameters
 848 
c     Auxiliary variables and parameters
792 
      real len1,len2,len3
 849 
      real len1,len2,len3
793 
      real val1,val2,val3
 850 
      real val1,val2,val3
794 
      real vx1,vy1
 851 
      real vx1,vy1
795 
      real vx2,vy2
 852 
      real vx2,vy2
796 
      real pi
 853 
      real pi
797 
      parameter (pi=3.14159265359)
 854 
      parameter (pi=3.14159265359)
798 
 
 855 
 
799 
      vx1 = ux1
 856 
      vx1 = ux1
800 
      vx2 = ux2
 857 
      vx2 = ux2
801 
      vy1 = uy1
 858 
      vy1 = uy1
802 
      vy2 = uy2
 859 
      vy2 = uy2
803 
 
 860 
 
804 
      len1=sqrt(vx1*vx1+vy1*vy1)
 861 
      len1=sqrt(vx1*vx1+vy1*vy1)
805 
      len2=sqrt(vx2*vx2+vy2*vy2)
 862 
      len2=sqrt(vx2*vx2+vy2*vy2)
806 
 
 863 
 
807 
      if ((len1.gt.0.).and.(len2.gt.0.)) then
 864 
      if ((len1.gt.0.).and.(len2.gt.0.)) then
808 
         vx1=vx1/len1
 865 
         vx1=vx1/len1
809 
         vy1=vy1/len1
 866 
         vy1=vy1/len1
810 
         vx2=vx2/len2
 867 
         vx2=vx2/len2
811 
         vy2=vy2/len2
 868 
         vy2=vy2/len2
812 
 
 869 
 
813 
         val1=vx1*vx2+vy1*vy2
 870 
         val1=vx1*vx2+vy1*vy2
814 
         val2=-vy1*vx2+vx1*vy2
 871 
         val2=-vy1*vx2+vx1*vy2
815 
 
 872 
 
816 
         len3=sqrt(val1*val1+val2*val2)
 873 
         len3=sqrt(val1*val1+val2*val2)
817 
 
 874 
 
818 
         if ( (val1.ge.0.).and.(val2.ge.0.) ) then
 875 
         if ( (val1.ge.0.).and.(val2.ge.0.) ) then
819 
            val3=acos(val1/len3)
 876 
            val3=acos(val1/len3)
820 
         else if ( (val1.lt.0.).and.(val2.ge.0.) ) then
 877 
         else if ( (val1.lt.0.).and.(val2.ge.0.) ) then
821 
            val3=pi-acos(abs(val1)/len3)
 878 
            val3=pi-acos(abs(val1)/len3)
822 
         else if ( (val1.ge.0.).and.(val2.le.0.) ) then
 879 
         else if ( (val1.ge.0.).and.(val2.le.0.) ) then
823 
            val3=-acos(val1/len3)
 880 
            val3=-acos(val1/len3)
824 
         else if ( (val1.lt.0.).and.(val2.le.0.) ) then
 881 
         else if ( (val1.lt.0.).and.(val2.le.0.) ) then
825 
            val3=-pi+acos(abs(val1)/len3)
 882 
            val3=-pi+acos(abs(val1)/len3)
826 
         endif
 883 
         endif
827 
      else
 884 
      else
828 
         val3=0.
 885 
         val3=0.
829 
      endif
 886 
      endif
830 
 
 887 
 
831 
      angle=180./pi*val3
 888 
      angle=180./pi*val3
832 
 
 889 
 
833 
      END
 890 
      END
834 
 
 891 
 
835 
 
 892 
 
836 
c     --------------------------------------------------------------------------
 893 
c     --------------------------------------------------------------------------
837 
c     Spherical distance between lat/lon points
894 
c     Spherical distance between lat/lon points
838 
c     --------------------------------------------------------------------------
 895 
c     --------------------------------------------------------------------------
839 
 
 896 
 
840 
      real function sdis(xp,yp,xq,yq)
 897 
      real function sdis(xp,yp,xq,yq)
841 
c
 898 
c
842 
c     calculates spherical distance (in km) between two points given
 899 
c     calculates spherical distance (in km) between two points given
843 
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
 900 
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
844 
c
 901 
c
845 
      real      re
 902 
      real      re
846 
      parameter (re=6370.)
 903 
      parameter (re=6370.)
847 
      real      pi180
 904 
      real      pi180
848 
      parameter (pi180=3.14159/180.)
 905 
      parameter (pi180=3.14159/180.)
849 
      real      xp,yp,xq,yq,arg
 906 
      real      xp,yp,xq,yq,arg
850 
 
 907 
 
851 
      arg=sin(pi180*yp)*sin(pi180*yq)+
 908 
      arg=sin(pi180*yp)*sin(pi180*yq)+
852 
     >    cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(xp-xq))
 909 
     >    cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(xp-xq))
853 
      if (arg.lt.-1.) arg=-1.
 910 
      if (arg.lt.-1.) arg=-1.
854 
      if (arg.gt.1.) arg=1.
 911 
      if (arg.gt.1.) arg=1.
855 
 
 912 
 
856 
      sdis=re*acos(arg)
 913 
      sdis=re*acos(arg)
857 
 
 914 
 
858 
      end
 915 
      end
859 
 
 916