#include "dims.h"
      SUBROUTINE SWDOT
      use cons_module,only: kut,kmax,imax,imaxp2,len1,kmaxp1,len3,
     |  expz,expzmid,boltz,p0,rmassinv_o2,rmassinv_o,rmassinv_n2
      implicit none
C     ****     CALCULATES W ON LINE NJ AND Z TO NJNP
#include "params.h"
#include "fcom.h"
#include "vscr.h"
#include "index.h"
#include "buff.h"
#include "phys.h"
!
! Local:
      integer :: i,k,kk
      integer :: ntk,nwk,nznpk,nzk,nps1k,nps2k,nnmbark,nmshk
      real :: fmin,fmax
      real :: wn(zimxp,zkmxp),pzps(zimxp,zkmxp)
!
C     ****     S1=DIVRG(V)
      CALL DIVRG(S1)
C     ****      INSERT PERIODIC POINTS  FOR W
      DO 4 I=1,2
      DO 4 K=1,KMAX
      S1(I,K)=S1(I+IMAX,K)
      S1(I+IMAXP2,K)=S1(I+2,K)
    4 CONTINUE
C     ****     W(KMAXP1)=S1(KMAX)
      NWK=NJNP+NW
      DO 1 I=1,LEN1
      F(I,NWK+KMAX)=S1(I,KMAX)
    1 CONTINUE
C     ****     W(K)=expzmid*(expzmid*W(K+1)+dz*S1(K))
      NWK=NWK+KMAX
      DO 22 KK=1,KMAX
      K=KMAXP1-KK
      NWK=NWK-1
      DO 22 I=1,LEN1
      F(I,NWK)=expzmid*(expzmid*F(I,NWK+1)+dz*S1(I,K))
   22 CONTINUE
      NZNPK=NJNP+NZ
      NZK=NJ+NZ
      DO 3 I=1,LEN3
      F(I,NZNPK)=F(I,NZK)
    3 CONTINUE
      NWK=NJNP+NW
      CALL FILTER(NWK,KMAXP1,KUT(J))
!
! Save vertical velocity:
      do i=1,len3
        wn(i,1) = 
      enddo  
C     ****
C     ****     Calculate n*M = p0*exp(-s)*M/(k*T)          (k+1/2)
C     ****      and  M = 1/(psi1/rm1+psi2/rm2+(1-psi1-psi2)/rm3)
C     ****      where:                                     (k+1/2)
C     ****
C     ****      M = mean molecular mass
C     ****      n = total number density
C     ****      p0 = standard pressure
C     ****      s = vertical pressure coordinate
C     ****      k = Boltzmann constant
C     ****      T = temperature
C     ****
      NPS1K = NJ+NPS-1
      NPS2K = NJ+NPS2-1
      NTK = NJ+NT-1
      NNMBARK = NNMBAR-1
      NMSHK = NMSH-1
      DO 5 K = 1,KMAX
        NPS1K = NPS1K+1
        NPS2K = NPS2K+1
        NTK = NTK+1
        NNMBARK = NNMBARK+1
        NMSHK = NMSHK+1
        DO I = 1,LEN1
          F(I,NMSHK) = 1./(F(I,NPS1K)*rmassinv_o2+F(I,NPS2K)*
     |      rmassinv_o+(1.-F(I,NPS1K)-F(I,NPS2K))*rmassinv_n2)
          F(I,NNMBARK) = p0*expz(K)*F(I,NMSHK)/(boltz*F(I,NTK))
        enddo
    5 CONTINUE
!     call fminmax(f(1,njnp+nz),zimxp*zkmxp,fmin,fmax)
!     write(6,"('swdot j=',i2,' f(1,njnp+nz) min,max=',2e12.4)") 
!    |  j,fmin,fmax
!     call fminmax(f(1,njnp+nw),zimxp*zkmxp,fmin,fmax)
!     write(6,"('swdot j=',i2,' f(1,njnp+nw) min,max=',2e12.4)") 
!    |  j,fmin,fmax
!     call fminmax(f(1,nnmbar),zimxp*zkmxp,fmin,fmax)
!     write(6,"('swdot j=',i2,' f(1,nnmbar) min,max=',2e12.4)") 
!    |  j,fmin,fmax
!     call fminmax(f(1,nmsh),zimxp*zkmxp,fmin,fmax)
!     write(6,"('swdot j=',i2,' f(1,nmsh) min,max=',2e12.4)") 
!    |  j,fmin,fmax
      RETURN
      END