program KINEMATIC c c this program computes the kinematic vertical motion using a tria c -ngular method . the observations need to be stratified in the ve c -rtical each layer must have at least three observations not on a c straight line . the layers are seperated by a check of alat=999. c the subroutine has 10 levels and a maximum of 75 point values of c direction and speed measured in degrees and knots ,respectively. c these observations are irregularly spaced . c the vertical motion, w is computed for the 10 levels using the c mean divergence of the layer.omega at the surface is set to zero c whereas at the top it is computed assuming that the layer diverg c -ence is the same as the layer below it. c c definitions of parameters c c nlvl : number of stratified layers c nrec : number of obs in the vertical layer (max 75) c pres : pressure level in mb (f6.1) c idir : wind direction in degrees (i3) c ispd : wind speed (knots,i3) c xo : longitude of the origin in degrees c yo : latitude of the origin in degrees c alat : latitudinal position of obs (f5.1) c alon : longitudinal position of obs (f5.1) c au,bu,cu : coefficients of planar surface of x-component c av,bv,cv : coefficients of planar surface of y-component c div : mean layer divergence (per sec) c vort : mean layer vorticity (per sec) c w : vertical velocity (mb/sec) c fact : conversion factor from knots to meter per sec c rpd : conversion factor from degrees to radians c npt : number of observations in a layer c divs : vertical integral of the mean divergence c diva : magnitude of the mean divergence c eps : correction factor of the mean divergence c dp : depth of the pressure levels (here constant) c parameter (nlvl = 8 ,nrec = 75) c c declaration of varibles c real div(8), vort(8), w(10), p1(10), ucomp(75) , vcomp(75) real ylat(75),xlong(75),au(8),bu(8),cu(8),av(8),bv(8),cv(8) c open(5,file='kinematic.dat',status='old') c data fact,xo,yo,dp/0.514791,42.,15.,100./ pi = 4.0*atan(1.0) divj = 0. divi = 0. rpd = pi/180. c c initialize the vertical motion to zero c do 3100 nk = 1, 10 3100 w(nk) = 0.0 c c loop over levels c do 3102 nk = 1, nlvl npt = 0 do 3104 n = 1, nrec c c read in observations with specified format c read(5,15,end=22)alat,alon,pres,idir,ispd 24 continue if (alat.eq.999.) go to 13 spd = float(ispd)*fact beta = float(idir)*rpd ucomp(n) = -spd*sin(beta) vcomp(n) = -spd*cos(beta) ylat (n) = (alat-yo)*1.111e5 xlong(n) = (alon-xo)*1.111e5*cos(alat)*rpd npt = npt+1 3104 continue 13 continue c c subroutine PLNSFC is called for both u- and v-component c to solve the linear system eq(3.3 to 3.5) in each layer . c call PLNSFC (ucomp,xlong,ylat,au(nk),bu(nk),cu(nk),npt) call PLNSFC (vcomp,xlong,ylat,av(nk),bv(nk),cv(nk),npt) 3102 continue 22 continue c c compute divergence and vorticity in the layer vorticity is added c just for complement .eqs(3.6 and 3.7) are used . c do 3106 k = 1, nlvl vort(k)= av(k)-bu(k) div(k) = au(k)+bv(k) divs = divs+div(k) diva = diva+abs(div(k)) 3106 continue c c compute the divergence in the ninth layer . c divs = divs+div(nlvl) diva = diva+abs(div(nlvl)) c c correct the divergence of the layer . c eps = -(divs/diva) do 3108 k = 1, nlvl div(k) = div(k)+eps*abs(div(k)) 3108 continue c c compute the omegas at all levels c w(1) = 0. do 3110 k = 2, 9 3110 w(k) = w(k-1)+div(k-1)*dp w(nlvl+2) = w(nlvl+1)+div(nlvl)*dp c c write output for the 8 first levels. c the format shows that the divergence c and vorticity are computed in the layers c between two pressure levels. c do 3112 k = 1, 10 3112 p1(k) = 1000.-((k-1)*100.) write(6,1000) write(6,1001) do 3114 k = 1, nlvl write (6,1002) p1(k) , w(k) if(k.ge.8) goto 3114 write(6,1003) div(k), vort(k) 3114 continue 15 format(2f5.1,f6.1,1x,2i3) 1000 format(//,1x,'pres lvl',4x,' omega ',6x,'divergence', + 5x,'vorticity') 1001 format(3x,'(mb)',5x,'(mb/sec)',6x,'(per sec)', + 6x,'(per sec)',/) 1002 format(1x,/,2x,f5.0,2x,e12.4) 1003 format(24x,e12.4,2x,e12.4) stop end subroutine PLNSFC (q,x,y,a,b,c,n) c c This subroutine solves the regression equations c by the least square approach to determine the c coefficients . t involves finding the solution c of the normal equations c real q(75),x(75),y(75),sum(8) do 3110 k = 1, 8 3110 sum(k) = 0 . do 3112 k = 1,n sum(1) = sum(1)+x(k) sum(2) = sum(2)+y(k) sum(3) = sum(3)+x(k)*x(k) sum(4) = sum(4)+x(k)*y(k) sum(5) = sum(5)+y(k)*y(k) sum(6) = sum(6)+q(k) sum(7) = sum(7)+q(k)*x(k) sum(8) = sum(8)+q(k)*y(k) 3112 continue denom = n*(sum(3)*sum(5)-sum(4)*sum(4)) & -(sum(1)*(sum(1)*sum(5)-sum(2)*sum(4))) & +(sum(2)*(sum(1)*sum(4)-sum(2)*sum(3))) c c = (sum(6)*(sum(3)*sum(5)-sum(4)*sum(4))- & (sum(7)*(sum(1)*sum(5)-sum(2)*sum(4)))+ & (sum(8)*(sum(1)*sum(4)-sum(2)*sum(3))))/denom c a = (n*(sum(7)*sum(5)-sum(4)*sum(8))- & (sum(1)*(sum(6)*sum(5)-sum(2)*sum(8)))+ & (sum(2)*(sum(6)*sum(4)-sum(2)*sum(7))))/denom c b = (n*(sum(3)*sum(8)-sum(7)*sum(4))- & ((sum(1)*(sum(1)*sum(8)-sum(6)*sum(4)))+ & (sum(2)*sum(1)*sum(7)-sum(6)*sum(3))))/denom c return end