Subversion Repositories lagranto.um

Rev

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

Rev 3 Rev 12
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,lon,
 554 
      subroutine calc_ABSVORT (absvort,dudy,dvdx,u,lat,lon,
555 
     >                                 pollon,pollat)
 555 
     >                                 pollon,pollat)
556 
 
 556 
 
557 
      implicit none
 557 
      implicit none
558 
 
 558 
 
559 
c     Argument declaration
 559 
c     Argument declaration
560 
      real  absvort       ! Absolute vorticity [s^-1]
 560 
      real  absvort       ! Absolute vorticity [s^-1]
561 
      real  u             ! Zonal wind component [m/s]
 561 
      real  u             ! Zonal wind component [m/s]
562 
      real  dudy          ! du/dy [s^-1]
 562 
      real  dudy          ! du/dy [s^-1]
563 
      real  dvdx          ! dv/dx [s^-1]
 563 
      real  dvdx          ! dv/dx [s^-1]
564 
      real  lat           ! Latitude [deg]
 564 
      real  lat           ! Latitude [deg]
565 
      real  lon           ! Longitude [deg]
 565 
      real  lon           ! Longitude [deg]
566 
      real  pollon,pollat ! Rotated north pole [deg]
 566 
      real  pollon,pollat ! Rotated north pole [deg]
567 
 
 567 
 
568 
c     Physical parameters
 568 
c     Physical parameters
569 
      real      pi180
 569 
      real      pi180
570 
      parameter (pi180=3.14159/180.)
 570 
      parameter (pi180=3.14159/180.)
571 
      real      deltay
 571 
      real      deltay
572 
      parameter (deltay =1.11e5)
 572 
      parameter (deltay =1.11e5)
573 
      real      omega
 573 
      real      omega
574 
      parameter (omega=7.292e-5)
 574 
      parameter (omega=7.292e-5)
575 
 
 575 
 
576 
c     Externals
576 
c     Externals
577 
      real      lmtolms
577 
      real      lmtolms
578 
      real      phtophs
578 
      real      phtophs
579 
      real      lmstolm
 579 
      real      lmstolm
580 
      real      phstoph
580 
      real      phstoph
581 
      external  lmtolms,phtophs
 581 
      external  lmtolms,phtophs
582 
      external  lmstolm,phstoph
 582 
      external  lmstolm,phstoph
583 
 
 583 
 
584 
c     Auxiliary varaibles
 584 
c     Auxiliary varaibles
585 
      real      truelat
 585 
      real      truelat
586 
 
 586 
 
587 
      truelat = phstoph(lat,lon,pollat,pollon)
 587 
      truelat = phstoph(lat,lon,pollat,pollon)
588 
 
 588 
 
589 
      absvort = dvdx - dudy + u * pi180/deltay * tan(pi180 * lat) +
 589 
      absvort = dvdx - dudy + u * pi180/deltay * tan(pi180 * lat) +
590 
     >          2. * omega * sin(pi180 * truelat)
590 
     >          2. * omega * sin(pi180 * truelat)
591
591
592 
      end
 592 
      end
593 
 
 593 
 
594 
c     -------------------------------------------------------------------------
 594 
c     -------------------------------------------------------------------------
595 
c     Divergence (DIV)
 595 
c     Divergence (DIV)
596 
c     -------------------------------------------------------------------------
 596 
c     -------------------------------------------------------------------------
597 
 
 597 
 
598 
      subroutine calc_DIV (div,dudx,dvdy,v,lat)
 598 
      subroutine calc_DIV (div,dudx,dvdy,v,lat)
599 
 
 599 
 
600 
      implicit none
 600 
      implicit none
601
601
602 
c     Argument declaration
 602 
c     Argument declaration
603 
      real  div        ! Divergence [s^-1]
 603 
      real  div        ! Divergence [s^-1]
604 
      real  v          ! Meridional wind component [m/s]
 604 
      real  v          ! Meridional wind component [m/s]
605 
      real  dudx       ! du/dx [s^-1]
 605 
      real  dudx       ! du/dx [s^-1]
606 
      real  dvdy       ! dv/dy [s^-1]
 606 
      real  dvdy       ! dv/dy [s^-1]
607 
      real  lat        ! Latitude [deg]
 607 
      real  lat        ! Latitude [deg]
608 
 
 608 
 
609 
c     Physical parameters
 609 
c     Physical parameters
610 
      real      pi180
 610 
      real      pi180
611 
      parameter (pi180=3.14159/180.)
 611 
      parameter (pi180=3.14159/180.)
612 
      real      deltay
 612 
      real      deltay
613 
      parameter (deltay =1.11e5)
 613 
      parameter (deltay =1.11e5)
614 
 
 614 
 
615 
      div = dudx + dvdy - v * pi180/deltay * tan(pi180 * lat)
615 
      div = dudx + dvdy - v * pi180/deltay * tan(pi180 * lat)
616
616
617 
      end
 617 
      end
618 
 
 618 
 
619 
c     -------------------------------------------------------------------------
 619 
c     -------------------------------------------------------------------------
620 
c     Deformation (DEF)
 620 
c     Deformation (DEF)
621 
c     -------------------------------------------------------------------------
 621 
c     -------------------------------------------------------------------------
622 
 
 622 
 
623 
      subroutine calc_DEF (def,dudx,dvdx,dudy,dvdy)
 623 
      subroutine calc_DEF (def,dudx,dvdx,dudy,dvdy)
624 
 
 624 
 
625 
      implicit none
 625 
      implicit none
626 
 
 626 
 
627 
c     Argument declaration
 627 
c     Argument declaration
628 
      real  def        ! Deformation [s^-1]
 628 
      real  def        ! Deformation [s^-1]
629 
      real  dudx       ! du/dx [s^-1]
 629 
      real  dudx       ! du/dx [s^-1]
630 
      real  dvdx       ! dv/dy [s^-1]
 630 
      real  dvdx       ! dv/dy [s^-1]
631 
      real  dudy       ! du/dx [s^-1]
 631 
      real  dudy       ! du/dx [s^-1]
632 
      real  dvdy       ! dv/dy [s^-1]
 632 
      real  dvdy       ! dv/dy [s^-1]
633 
 
 633 
 
634 
c     Physical parameters
 634 
c     Physical parameters
635 
      real      pi180
 635 
      real      pi180
636 
      parameter (pi180=3.14159/180.)
 636 
      parameter (pi180=3.14159/180.)
637 
 
 637 
 
638 
      def = sqrt( (dvdx+dudy)**2 + (dudx-dvdy)**2 )
 638 
      def = sqrt( (dvdx+dudy)**2 + (dudx-dvdy)**2 )
639
639
640 
      end
 640 
      end
641 
 
 641 
 
642 
c     -------------------------------------------------------------------------
 642 
c     -------------------------------------------------------------------------
643 
c     Potential Vorticity (PV)
 643 
c     Potential Vorticity (PV)
644 
c     -------------------------------------------------------------------------
 644 
c     -------------------------------------------------------------------------
645 
 
 645 
 
646 
      subroutine calc_PV (pv,absvort,dthdp,dudp,dvdp,dthdx,dthdy)
 646 
      subroutine calc_PV (pv,absvort,dthdp,dudp,dvdp,dthdx,dthdy)
647 
 
 647 
 
648 
      implicit none
 648 
      implicit none
649 
 
 649 
 
650 
c     Argument declaration
 650 
c     Argument declaration
651 
      real  pv         ! Ertel-PV [PVU]
 651 
      real  pv         ! Ertel-PV [PVU]
652 
      real  absvort    ! Absolute vorticity [s^-1]
 652 
      real  absvort    ! Absolute vorticity [s^-1]
653 
      real  dthdp      ! dth/dp [K/Pa]
 653 
      real  dthdp      ! dth/dp [K/Pa]
654 
      real  dudp       ! du/dp [m/s per Pa]
 654 
      real  dudp       ! du/dp [m/s per Pa]
655 
      real  dvdp       ! dv/dp [m/s per Pa]
 655 
      real  dvdp       ! dv/dp [m/s per Pa]
656 
      real  dthdx      ! dth/dx [K/m]
 656 
      real  dthdx      ! dth/dx [K/m]
657 
      real  dthdy      ! dth/dy [K/m]
 657 
      real  dthdy      ! dth/dy [K/m]
658 
 
 658 
 
659 
c     Physical and numerical parameters
 659 
c     Physical and numerical parameters
660 
      real      scale
 660 
      real      scale
661 
      parameter (scale=1.E6)
 661 
      parameter (scale=1.E6)
662 
      real      g
 662 
      real      g
663 
      parameter (g=9.80665)
 663 
      parameter (g=9.80665)
664 
 
 664 
 
665 
      pv = -scale * g * ( absvort * dthdp + dudp * dthdy - dvdp * dthdx)
 665 
      pv = -scale * g * ( absvort * dthdp + dudp * dthdy - dvdp * dthdx)
666 
 
 666 
 
667 
      end
 667 
      end
668 
 
 668 
 
669 
c     -------------------------------------------------------------------------
 669 
c     -------------------------------------------------------------------------
670 
c     Richardson number (RI)
 670 
c     Richardson number (RI)
671 
c     -------------------------------------------------------------------------
 671 
c     -------------------------------------------------------------------------
672 
 
 672 
 
673 
      subroutine calc_RI (ri,dudp,dvdp,nsq,rho)
 673 
      subroutine calc_RI (ri,dudp,dvdp,nsq,rho)
674 
 
 674 
 
675 
      implicit none
 675 
      implicit none
676 
 
 676 
 
677 
c     Argument declaration
 677 
c     Argument declaration
678 
      real  ri         ! Richardson number
 678 
      real  ri         ! Richardson number
679 
      real  dudp       ! Du/Dp [m/s per Pa]
 679 
      real  dudp       ! Du/Dp [m/s per Pa]
680 
      real  dvdp       ! Dv/Dp [m/s per Pa]
 680 
      real  dvdp       ! Dv/Dp [m/s per Pa]
681 
      real  nsq        ! Squared Brunt-Vailälä frequency [s^-1]
 681 
      real  nsq        ! Squared Brunt-Vailälä frequency [s^-1]
682 
      real  rho        ! Density [kg/m^3]
 682 
      real  rho        ! Density [kg/m^3]
683 
 
 683 
 
684 
c     Physical and numerical parameters
 684 
c     Physical and numerical parameters
685 
      real      g
 685 
      real      g
686 
      parameter (g=9.80665)
 686 
      parameter (g=9.80665)
687 
 
 687 
 
688 
      ri = nsq / ( dudp**2 + dvdp**2 ) / ( rho * g )**2
 688 
      ri = nsq / ( dudp**2 + dvdp**2 ) / ( rho * g )**2
689 
 
 689 
 
690 
      end
 690 
      end
691 
 
 691 
 
692 
c     -------------------------------------------------------------------------
 692 
c     -------------------------------------------------------------------------
693 
c     Ellrod and Knapp's turbulence indicator (TI)
 693 
c     Ellrod and Knapp's turbulence indicator (TI)
694 
c     -------------------------------------------------------------------------
 694 
c     -------------------------------------------------------------------------
695 
 
 695 
 
696 
      subroutine calc_TI (ti,def,dudp,dvdp,rho)
 696 
      subroutine calc_TI (ti,def,dudp,dvdp,rho)
697
697
698 
      implicit none
 698 
      implicit none
699 
 
 699 
 
700 
c     Argument declaration
 700 
c     Argument declaration
701 
      real  ti         ! Turbulence idicator
 701 
      real  ti         ! Turbulence idicator
702 
      real  def        ! Deformation [s^-1]
 702 
      real  def        ! Deformation [s^-1]
703 
      real  dudp       ! Du/Dp [m/s per Pa]
 703 
      real  dudp       ! Du/Dp [m/s per Pa]
704 
      real  dvdp       ! Dv/Dp [m/s per Pa]
 704 
      real  dvdp       ! Dv/Dp [m/s per Pa]
705 
      real  rho        ! Density [kg/m^3]
 705 
      real  rho        ! Density [kg/m^3]
706 
 
 706 
 
707 
c     Physical and numerical parameters
 707 
c     Physical and numerical parameters
708 
      real      g
 708 
      real      g
709 
      parameter (g=9.80665)
 709 
      parameter (g=9.80665)
710 
 
 710 
 
711 
      ti = def * sqrt ( dudp**2 + dvdp**2 ) * ( rho * g )
 711 
      ti = def * sqrt ( dudp**2 + dvdp**2 ) * ( rho * g )
712 
 
 712 
 
713 
      end
 713 
      end
714 
 
 714 
 
715 
c     -------------------------------------------------------------------------
 715 
c     -------------------------------------------------------------------------
716 
c     Distance from starting position
 716 
c     Distance from starting position
717 
c     -------------------------------------------------------------------------
 717 
c     -------------------------------------------------------------------------
718 
 
 718 
 
719 
      subroutine calc_DIST0 (dist0,lon0,lat0,lon1,lat1)
 719 
      subroutine calc_DIST0 (dist0,lon0,lat0,lon1,lat1)
720 
 
 720 
 
721 
      implicit none
 721 
      implicit none
722 
 
 722 
 
723 
c     Argument declaration
 723 
c     Argument declaration
724 
      real  dist0      ! Distance from starting position [km]
 724 
      real  dist0      ! Distance from starting position [km]
725 
      real  lon0,lat0  ! Starting position
 725 
      real  lon0,lat0  ! Starting position
726 
      real  lon1,lat1  ! New position
726 
      real  lon1,lat1  ! New position
727 
 
 727 
 
728 
c     Externals
728 
c     Externals
729 
      real     sdis
 729 
      real     sdis
730 
      external sdis
 730 
      external sdis
731 
 
 731 
 
732 
      dist0 = sdis(lon0,lat0,lon1,lat1)
 732 
      dist0 = sdis(lon0,lat0,lon1,lat1)
733
733
734 
      end
 734 
      end
735 
 
 735 
 
736 
c     -------------------------------------------------------------------------
 736 
c     -------------------------------------------------------------------------
737 
c     Heading of the trajectory (HEAD)
 737 
c     Heading of the trajectory (HEAD)
738 
c     -------------------------------------------------------------------------
 738 
c     -------------------------------------------------------------------------
739 
 
 739 
 
740 
      subroutine calc_HEAD (head,lon0,lat0,lon1,lat1)
 740 
      subroutine calc_HEAD (head,lon0,lat0,lon1,lat1)
741
741
742 
      implicit none
 742 
      implicit none
743 
 
 743 
 
744 
c     Argument declaration
 744 
c     Argument declaration
745 
      real  head       ! Heading angle (in deg) relativ to zonal direction
 745 
      real  head       ! Heading angle (in deg) relativ to zonal direction
746 
      real  lon0,lat0  ! Starting position
 746 
      real  lon0,lat0  ! Starting position
747 
      real  lon1,lat1  ! New position
747 
      real  lon1,lat1  ! New position
748 
 
 748 
 
749 
c     Physical parameters
 749 
c     Physical parameters
750 
      real      pi180
 750 
      real      pi180
751 
      parameter (pi180=3.14159/180.)
 751 
      parameter (pi180=3.14159/180.)
752 
 
 752 
 
753 
c     Auixiliary variables
 753 
c     Auixiliary variables
754 
      real     dx,dy
 754 
      real     dx,dy
755 
 
 755 
 
756 
      dx = (lon1-lon0) * cos(pi180*0.5*(lat0+lat1))
 756 
      dx = (lon1-lon0) * cos(pi180*0.5*(lat0+lat1))
757 
      dy = lat1-lat0
 757 
      dy = lat1-lat0
758 
 
 758 
 
759 
      call getangle(1.,0.,dx,dy,head)
 759 
      call getangle(1.,0.,dx,dy,head)
760 
 
 760 
 
761 
      end
 761 
      end
762 
 
 762 
 
763 
c
763 
c
764 
c     *************************************************************************
 764 
c     *************************************************************************
765 
c     Auxiliary subroutines and functions
 765 
c     Auxiliary subroutines and functions
766 
c     *************************************************************************
 766 
c     *************************************************************************
767
767
768 
c     -------------------------------------------------------------------------
 768 
c     -------------------------------------------------------------------------
769 
c     Saturation vapor pressure over water
 769 
c     Saturation vapor pressure over water
770 
c     -------------------------------------------------------------------------
 770 
c     -------------------------------------------------------------------------
771 
 
 771 
 
772 
      real function esat(t)
 772 
      real function esat(t)
773
773
774 
C     This function returns the saturation vapor pressure over water
 774 
C     This function returns the saturation vapor pressure over water
775 
c     (mb) given the temperature (Kelvin).
 775 
c     (mb) given the temperature (Kelvin).
776 
C     The algorithm is due to Nordquist, W. S. ,1973: "Numerical
 776 
C     The algorithm is due to Nordquist, W. S. ,1973: "Numerical
777 
C     Approximations of Selected Meteorological Parameters for Cloud
 777 
C     Approximations of Selected Meteorological Parameters for Cloud
778 
C     Physics Problems" ECOM-5475, Atmospheric Sciences Laboratory,
 778 
C     Physics Problems" ECOM-5475, Atmospheric Sciences Laboratory,
779 
c     U. S. Army Electronics Command, White Sands Missile Range,
 779 
c     U. S. Army Electronics Command, White Sands Missile Range,
780 
c     New Mexico 88002.
 780 
c     New Mexico 88002.
781
781
782 
      real p1,p2,c1,t
 782 
      real p1,p2,c1,t
783
783
784 
      p1=11.344-0.0303998*t
 784 
      p1=11.344-0.0303998*t
785 
      p2=3.49149-1302.8844/t
 785 
      p2=3.49149-1302.8844/t
786 
      c1=23.832241-5.02808*log10(t)
 786 
      c1=23.832241-5.02808*log10(t)
787 
      esat=10.**(c1-1.3816e-7*10.**p1+8.1328e-3*10.**p2-2949.076/t)
 787 
      esat=10.**(c1-1.3816e-7*10.**p1+8.1328e-3*10.**p2-2949.076/t)
788 
 
 788 
 
789 
      end
 789 
      end
790 
 
 790 
 
791 
c     --------------------------------------------------------------------------
 791 
c     --------------------------------------------------------------------------
792 
c     Angle between two vectors
 792 
c     Angle between two vectors
793 
c     --------------------------------------------------------------------------
 793 
c     --------------------------------------------------------------------------
794 
 
 794 
 
795 
      SUBROUTINE getangle (ux1,uy1,ux2,uy2,angle)
 795 
      SUBROUTINE getangle (ux1,uy1,ux2,uy2,angle)
796 
 
 796 
 
797 
c     Given two vectors <ux1,uy1> and <ux2,uy2>, determine the angle (in deg)
 797 
c     Given two vectors <ux1,uy1> and <ux2,uy2>, determine the angle (in deg)
798 
c     between the two vectors.
 798 
c     between the two vectors.
799 
 
 799 
 
800 
      implicit none
 800 
      implicit none
801 
 
 801 
 
802 
c     Declaration of subroutine parameters
 802 
c     Declaration of subroutine parameters
803 
      real ux1,uy1
 803 
      real ux1,uy1
804 
      real ux2,uy2
 804 
      real ux2,uy2
805 
      real angle
 805 
      real angle
806 
 
 806 
 
807 
c     Auxiliary variables and parameters
 807 
c     Auxiliary variables and parameters
808 
      real len1,len2,len3
 808 
      real len1,len2,len3
809 
      real val1,val2,val3
 809 
      real val1,val2,val3
810 
      real vx1,vy1
 810 
      real vx1,vy1
811 
      real vx2,vy2
 811 
      real vx2,vy2
812 
      real pi
 812 
      real pi
813 
      parameter (pi=3.14159265359)
 813 
      parameter (pi=3.14159265359)
814 
 
 814 
 
815 
      vx1 = ux1
 815 
      vx1 = ux1
816 
      vx2 = ux2
 816 
      vx2 = ux2
817 
      vy1 = uy1
 817 
      vy1 = uy1
818 
      vy2 = uy2
 818 
      vy2 = uy2
819 
 
 819 
 
820 
      len1=sqrt(vx1*vx1+vy1*vy1)
 820 
      len1=sqrt(vx1*vx1+vy1*vy1)
821 
      len2=sqrt(vx2*vx2+vy2*vy2)
 821 
      len2=sqrt(vx2*vx2+vy2*vy2)
822 
 
 822 
 
823 
      if ((len1.gt.0.).and.(len2.gt.0.)) then
 823 
      if ((len1.gt.0.).and.(len2.gt.0.)) then
824 
         vx1=vx1/len1
 824 
         vx1=vx1/len1
825 
         vy1=vy1/len1
 825 
         vy1=vy1/len1
826 
         vx2=vx2/len2
 826 
         vx2=vx2/len2
827 
         vy2=vy2/len2
 827 
         vy2=vy2/len2
828 
 
 828 
 
829 
         val1=vx1*vx2+vy1*vy2
 829 
         val1=vx1*vx2+vy1*vy2
830 
         val2=-vy1*vx2+vx1*vy2
 830 
         val2=-vy1*vx2+vx1*vy2
831 
 
 831 
 
832 
         len3=sqrt(val1*val1+val2*val2)
 832 
         len3=sqrt(val1*val1+val2*val2)
833 
 
 833 
 
834 
         if ( (val1.ge.0.).and.(val2.ge.0.) ) then
 834 
         if ( (val1.ge.0.).and.(val2.ge.0.) ) then
835 
            val3=acos(val1/len3)
 835 
            val3=acos(val1/len3)
836 
         else if ( (val1.lt.0.).and.(val2.ge.0.) ) then
 836 
         else if ( (val1.lt.0.).and.(val2.ge.0.) ) then
837 
            val3=pi-acos(abs(val1)/len3)
 837 
            val3=pi-acos(abs(val1)/len3)
838 
         else if ( (val1.ge.0.).and.(val2.le.0.) ) then
 838 
         else if ( (val1.ge.0.).and.(val2.le.0.) ) then
839 
            val3=-acos(val1/len3)
 839 
            val3=-acos(val1/len3)
840 
         else if ( (val1.lt.0.).and.(val2.le.0.) ) then
 840 
         else if ( (val1.lt.0.).and.(val2.le.0.) ) then
841 
            val3=-pi+acos(abs(val1)/len3)
 841 
            val3=-pi+acos(abs(val1)/len3)
842 
         endif
 842 
         endif
843 
      else
 843 
      else
844 
         val3=0.
 844 
         val3=0.
845 
      endif
 845 
      endif
846 
 
 846 
 
847 
      angle=180./pi*val3
 847 
      angle=180./pi*val3
848 
 
 848 
 
849 
      END
 849 
      END
850 
 
 850 
 
851 
 
 851 
 
852 
c     --------------------------------------------------------------------------
 852 
c     --------------------------------------------------------------------------
853 
c     Spherical distance between lat/lon points
853 
c     Spherical distance between lat/lon points
854 
c     --------------------------------------------------------------------------
 854 
c     --------------------------------------------------------------------------
855 
 
 855 
 
856 
      real function sdis(xp,yp,xq,yq)
 856 
      real function sdis(xp,yp,xq,yq)
857 
c
 857 
c
858 
c     calculates spherical distance (in km) between two points given
 858 
c     calculates spherical distance (in km) between two points given
859 
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
 859 
c     by their spherical coordinates (xp,yp) and (xq,yq), respectively.
860 
c
 860 
c
861 
      real      re
 861 
      real      re
862 
      parameter (re=6370.)
 862 
      parameter (re=6370.)
863 
      real      pi180
 863 
      real      pi180
864 
      parameter (pi180=3.14159/180.)
 864 
      parameter (pi180=3.14159/180.)
865 
      real      xp,yp,xq,yq,arg
 865 
      real      xp,yp,xq,yq,arg
866 
 
 866 
 
867 
      arg=sin(pi180*yp)*sin(pi180*yq)+
 867 
      arg=sin(pi180*yp)*sin(pi180*yq)+
868 
     >    cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(xp-xq))
 868 
     >    cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(xp-xq))
869 
      if (arg.lt.-1.) arg=-1.
 869 
      if (arg.lt.-1.) arg=-1.
870 
      if (arg.gt.1.) arg=1.
 870 
      if (arg.gt.1.) arg=1.
871 
 
 871 
 
872 
      sdis=re*acos(arg)
 872 
      sdis=re*acos(arg)
873 
 
 873 
 
874 
      end
 874 
      end
875 
 
 875