!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE GETLAPS                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE getlaps(nx_ext,ny_ext,nz_ext,dir_extd,                       & 1,8
           extdinit,extdfcst,julfname,i4time,                           &
           iproj_ext,scale_ext,                                         &
           trlon_ext,latnot_ext,x0_ext,y0_ext,                          &
           lat_ext,lon_ext,                                             &
           p_ext,hgt_ext,t_ext,qv_ext,u_ext,v_ext,                      &
           qc_ext,qr_ext,qi_ext,qs_ext,qh_ext,                          &
           istatus)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  OLAPS version.
!
!  Reads OLAPS file for processing by ext2arps, a program
!  that converts external files to ARPS variables and format.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Keith Brewster
!  May, 1994
!
!  MODIFICATION HISTORY:
!  September, 1994 (KB)
!  Upgrade for LLNL effort.
!
!  November, 1994 (KB)
!  Beefed up documentation.
!  OLAPS Version.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    dir_extd      Directory name for external file
!    extdinit      Initialized time in mm-dd-yyyy:hh:mm:ss format
!    extdfcst      Forecast hour in HHH:MM:SS format
!    julfname      File name in yyjjjhhmm format
!    i4time        Absolute time in seconds (from 1 Jan 1970) of desired data
!
!  OUTPUT:
!
!    iproj_ext     Map projection number of external data
!    scale_ext     Scale factor of external data
!    trlon_ext     True longitude of external data (degrees E)
!    latnot_ext(2) True latitude(s) of external data (degrees N)
!    x0_ext        x coordinate of origin of external data
!    y0_ext        y coordinate of origin of external data
!    lat_ext       latitude of external data points (degrees N)
!    lon_ext       longitude of external data points (degrees E)
!    p_ext         pressure (Pascal)
!    hgt_ext       height (m)
!    t_ext         temperature (K)
!    qv_ext        specific humidity (kg/kg)
!    u_ext         u wind component (m/s)
!    v_ext         v wind component (m/s)
!    qc_ext        Cloud water mixing ratio (kg/kg)
!    qr_ext        Rain water mixing ratio (kg/kg)
!    qi_ext        Ice mixing ratio (kg/kg)
!    qs_ext        Snow mixing ratio (kg/kg)
!    qh_ext        Hail mixing ratio (kg/kg)
!    istatus       status indicator
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  CHARACTER (LEN=*) :: dir_extd
  CHARACTER (LEN=19) :: extdinit
  CHARACTER (LEN=9) :: extdfcst
  CHARACTER (LEN=9) :: julfname
  INTEGER :: i4time
!
!-----------------------------------------------------------------------
!
!  Original grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj
  REAL :: scale,trlon,x0,y0
  REAL :: latnot(2)
!
!-----------------------------------------------------------------------
!
!  External grid variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: iproj_ext
  REAL :: scale_ext,trlon_ext
  REAL :: latnot_ext(2)
  REAL :: x0_ext,y0_ext
  REAL :: olatsw,olonsw
!
!-----------------------------------------------------------------------
!
!  Output external variable arrays
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx_ext,ny_ext,nz_ext

  REAL :: lat_ext(nx_ext,ny_ext)
  REAL :: lon_ext(nx_ext,ny_ext)
  REAL :: p_ext(nx_ext,ny_ext,nz_ext)     ! Pressure (Pascals)
  REAL :: hgt_ext(nx_ext,ny_ext,nz_ext)   ! Height (m)
  REAL :: t_ext(nx_ext,ny_ext,nz_ext)     ! Temperature (K)
  REAL :: qv_ext(nx_ext,ny_ext,nz_ext)    ! Specific humidity (kg/kg)
  REAL :: u_ext(nx_ext,ny_ext,nz_ext)     ! Eastward wind component
  REAL :: v_ext(nx_ext,ny_ext,nz_ext)     ! Northward wind component
  REAL :: qc_ext(nx_ext,ny_ext,nz_ext)    ! Cloud H2O mixing ratio (kg/kg)
  REAL :: qr_ext(nx_ext,ny_ext,nz_ext)    ! Rain  H2O mixing ratio (kg/kg)
  REAL :: qi_ext(nx_ext,ny_ext,nz_ext)    ! Ice   mixing ratio (kg/kg)
  REAL :: qs_ext(nx_ext,ny_ext,nz_ext)    ! Snow  mixing ratio (kg/kg)
  REAL :: qh_ext(nx_ext,ny_ext,nz_ext)    ! Hail  mixing ratio (kg/kg)

  INTEGER :: istatus
!
!-----------------------------------------------------------------------
!
!  OLAPS file grid variables
!
!-----------------------------------------------------------------------
!
  REAL :: dxlaps,dylaps
  PARAMETER (dxlaps=10000., dylaps=10000.)
!
  REAL :: xlaps(nx_ext),ylaps(ny_ext)
  REAL :: lat(nx_ext,ny_ext),lon(nx_ext,ny_ext)
  REAL :: terr(nx_ext,ny_ext)
!
!-----------------------------------------------------------------------
!
!  OLAPS Forecast variables
!
!-----------------------------------------------------------------------
!
  CHARACTER (LEN=50) :: dir_lsx,dir_lw3,dir_lq3,dir_lt1,dir_static
  CHARACTER (LEN=31) :: ext,ext_s

  INTEGER :: nz_ext_mx
  PARAMETER( nz_ext_mx = 41)
  INTEGER :: lvl(nz_ext_mx)
  CHARACTER (LEN=3) :: var_t(nz_ext_mx),var_q(nz_ext_mx)
  CHARACTER (LEN=3) :: var_u(nz_ext_mx),var_v(nz_ext_mx)
  CHARACTER (LEN=4) :: lvl_coord(nz_ext_mx)

  CHARACTER (LEN=4) :: sfc_coord
  CHARACTER (LEN=10) :: units(nz_ext)
  CHARACTER (LEN=125) :: comment(nz_ext)
!
  REAL :: sfct(nx_ext,ny_ext)
  REAL :: sfcmr(nx_ext,ny_ext)
  REAL :: sfcp_pa(nx_ext,ny_ext)
  REAL :: psnd(nz_ext)
  REAL :: tsnd(nz_ext)
  REAL :: htsnd(nz_ext)
  REAL :: mrsnd(nz_ext)
!
!-----------------------------------------------------------------------
!
!  OLAPS initializations for reading subroutines
!
!-----------------------------------------------------------------------
!
  DATA lvl /1100,1075,1050,1025,1000,975,950,925,900,                   &
                  875, 850, 825, 800,775,750,725,700,                   &
                  675, 650, 625, 600,575,550,525,500,                   &
                  475, 450, 425, 400,375,350,325,300,                   &
                  275, 250, 225, 200,175,150,125,100/
!
  DATA lvl_coord                                                        &
          /'HPA ','HPA ','HPA ','HPA ','HPA ','HPA ','HPA ','HPA ',     &
           'HPA ','HPA ','HPA ','HPA ','HPA ','HPA ','HPA ','HPA ',     &
           'HPA ','HPA ','HPA ','HPA ','HPA ','HPA ','HPA ','HPA ',     &
           'HPA ','HPA ','HPA ','HPA ','HPA ','HPA ','HPA ','HPA ',     &
           'HPA ','HPA ','HPA ','HPA ','HPA ','HPA ','HPA ','HPA ',     &
           'HPA '/
  DATA var_t /'t  ','t  ','t  ','t  ','t  ','t  ','t  ','t  ',          &
              't  ','t  ','t  ','t  ','t  ','t  ','t  ','t  ',          &
              't  ','t  ','t  ','t  ','t  ','t  ','t  ','t  ',          &
              't  ','t  ','t  ','t  ','t  ','t  ','t  ','t  ',          &
              't  ','t  ','t  ','t  ','t  ','t  ','t  ','t  ',          &
              't  '/
  DATA var_q /'sh ','sh ','sh ','sh ','sh ','sh ','sh ','sh ',          &
              'sh ','sh ','sh ','sh ','sh ','sh ','sh ','sh ',          &
              'sh ','sh ','sh ','sh ','sh ','sh ','sh ','sh ',          &
              'sh ','sh ','sh ','sh ','sh ','sh ','sh ','sh ',          &
              'sh ','sh ','sh ','sh ','sh ','sh ','sh ','sh ',          &
              'sh '/
  DATA var_u /'u  ','u  ','u  ','u  ','u  ','u  ','u  ','u  ',          &
              'u  ','u  ','u  ','u  ','u  ','u  ','u  ','u  ',          &
              'u  ','u  ','u  ','u  ','u  ','u  ','u  ','u  ',          &
              'u  ','u  ','u  ','u  ','u  ','u  ','u  ','u  ',          &
              'u  ','u  ','u  ','u  ','u  ','u  ','u  ','u  ',          &
              'u  '/
  DATA var_v /'v  ','v  ','v  ','v  ','v  ','v  ','v  ','v  ',          &
              'v  ','v  ','v  ','v  ','v  ','v  ','v  ','v  ',          &
              'v  ','v  ','v  ','v  ','v  ','v  ','v  ','v  ',          &
              'v  ','v  ','v  ','v  ','v  ','v  ','v  ','v  ',          &
              'v  ','v  ','v  ','v  ','v  ','v  ','v  ','v  ',          &
              'v  '/
!
!-----------------------------------------------------------------------
!
!  Misc internal variables
!
!-----------------------------------------------------------------------
!
  REAL :: smfact
  PARAMETER (smfact=0.5)
  INTEGER :: i,j,k,ldir,kdim,lsfc
  REAL :: swlat,swlon,xsw,ysw,grid_spacing
  REAL :: qvlast,t_sfc,w_sfc
!
!-----------------------------------------------------------------------
!
!  Functions
!
!-----------------------------------------------------------------------
!
  REAL :: qvtomxr
!
!-----------------------------------------------------------------------
!
!  Include files
!
!-----------------------------------------------------------------------
!
  INCLUDE 'lapsparms.cmn'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  iproj_ext=1          ! polar stereographic
  scale_ext=1.0        ! report lengths in m
  trlon_ext=-105.      ! orientation of OLAPS data grids
  latnot_ext(1)=40.
  latnot_ext(2)=40.
  x0_ext=0.
  y0_ext=0.
  swlat=32.5353
  swlon=-103.7612
!
  IF(dir_extd(1:5) == '     ') THEN
    dir_extd='/vortex/vortex95/lapsprd/'
    WRITE(6,'(a,a)')                                                    &
        ' Use the default directory of OLAPS analyses: ',               &
        ' /vortex/vortex95/lapsprd'
  ELSE
    WRITE(6,'(a,a)')                                                    &
        ' Use the input directory of OLAPS analyses: ',dir_extd
  END IF
!
  ldir=LEN(dir_extd)
  CALL strlnth(dir_extd,ldir)
  IF(dir_extd(ldir:ldir) /= '/') THEN
    ldir=ldir+1
    dir_extd(ldir:ldir)='/'
  END IF
!
  WRITE(dir_lsx,'(a,a)') dir_extd(1:ldir),'lsx/'
  WRITE(dir_lq3,'(a,a)') dir_extd(1:ldir),'lq3/'
  WRITE(dir_lw3,'(a,a)') dir_extd(1:ldir),'lw3/'
  WRITE(dir_lt1,'(a,a)') dir_extd(1:ldir),'lt1/'
!
  WRITE(6,'(a)')                                                        &
      ' The time in format yydddhhff of OLAPS initial file is ',julfname
!
  CALL cv_asc_i4time(julfname,i4time)
!
  DO k=1,nz_ext
    DO j=1,ny_ext
      DO i=1,nx_ext
        p_ext(i,j,k)=100.*FLOAT(lvl(k))
      END DO
    END DO
  END DO
!
  CALL getmapr(iproj,scale,latnot,trlon,x0,y0)
  CALL setmapr(iproj_ext,scale_ext,latnot_ext,trlon_ext)
  CALL setorig(1,x0_ext,y0_ext)
!
!-----------------------------------------------------------------------
!
!  Fake that get_laps_parms was called.
!  This is to insure proper reading of OLAPS terrain
!  file.
!
!-----------------------------------------------------------------------
!
  standard_latitude=latnot_ext(1)
  standard_longitude=trlon_ext
  iflag_lapsparms_cmn = 1
!
  dir_static='/vortex/vortex95/static/'
!
  ext_s='vortex95'
!
!-----------------------------------------------------------------------
!
!  Get OLAPS static variables
!
!-----------------------------------------------------------------------
!
  CALL rd_laps_static(dir_static,                                       &
               ext_s,nx_ext,ny_ext,1,                                   &
               'LAT',units,comment,                                     &
               lat_ext,grid_spacing,istatus)
  CALL rd_laps_static(dir_static,                                       &
               ext_s,nx_ext,ny_ext,1,                                   &
               'LON',units,comment,                                     &
               lon_ext,grid_spacing,istatus)
  CALL rd_laps_static(dir_static,                                       &
               ext_s,nx_ext,ny_ext,1,                                   &
               'AVG',units,comment,                                     &
               terr,grid_spacing,istatus)
!
!-----------------------------------------------------------------------
!
!  Reset map projection to previous values
!
!-----------------------------------------------------------------------
!
  CALL setmapr(iproj,scale,latnot,trlon)
  CALL setorig(1,x0,y0)
!
!-----------------------------------------------------------------------
!
!  Get file name information
!
!-----------------------------------------------------------------------
!
  ldir=LEN(dir_extd)
  CALL strlnth( dir_extd, ldir)
  IF(dir_extd(ldir:ldir) /= '/') THEN
    ldir=ldir+1
    dir_extd(ldir:ldir)='/'
  END IF
!
!-----------------------------------------------------------------------
!
!  Read LAPS surface pressure data
!
!-----------------------------------------------------------------------
!
  ext='lsx'
  kdim=1
  lsfc=0
  sfc_coord='AGL'
  CALL read_laps_data(i4time,dir_lsx,ext,                               &
                      nx_ext,ny_ext,1,kdim,                             &
                      'ps ',lsfc,sfc_coord,units,comment,               &
                      sfcp_pa,istatus)
  WRITE(6,'(a)')    ' Sfc Pressure Read'
  WRITE(6,'(a,i6)') ' istatus= ',istatus
  WRITE(6,'(a,a)')  ' units= ',units(1)
  IF(istatus /= 1) GO TO 598
!
!-----------------------------------------------------------------------
!
!  Read LAPS surface temperature data
!
!-----------------------------------------------------------------------
!
  CALL read_laps_data(i4time,dir_lsx,ext,nx_ext,ny_ext,1,kdim,          &
                      't  ',lsfc,sfc_coord,units,comment,               &
                      sfct,istatus)
  WRITE(6,'(a)')    ' Sfc Temp Read'
  WRITE(6,'(a,i6)') ' istatus= ',istatus
  IF(istatus /= 1) GO TO 598
!
!-----------------------------------------------------------------------
!
!  Read LAPS surface specific humidity
!
!-----------------------------------------------------------------------
!
  CALL read_laps_data(i4time,dir_lsx,ext,nx_ext,ny_ext,1,kdim,          &
                      'mr ',lsfc,sfc_coord,units,comment,               &
                      sfcmr,istatus)
  WRITE(6,'(a)')    ' Sfc specific humidity read'
  WRITE(6,'(a,i6)') ' istatus= ',istatus
  WRITE(6,'(a,a)')  ' units= ',units(1)
  IF(istatus /= 1) GO TO 598
!
!-----------------------------------------------------------------------
!
!  Read LAPS Temperature data
!
!-----------------------------------------------------------------------
!
  ext='lt1'
  kdim=nz_ext
  CALL read_laps_data(i4time,dir_lt1,ext,                               &
                      nx_ext,ny_ext,nz_ext,kdim,                        &
                      var_t,lvl,lvl_coord,units,comment,                &
                      t_ext, istatus)
  WRITE(6,'(a)')    ' Temperature Read'
  WRITE(6,'(a,i6)') ' istatus= ',istatus
  WRITE(6,'(a,a)')  ' units= ',units(1)
  IF(istatus /= 1) GO TO 598
!
!-----------------------------------------------------------------------
!
!  Read LAPS specific humidity field
!
!-----------------------------------------------------------------------
!
  ext='lq3'
  kdim=nz_ext
  CALL read_laps_data(i4time,dir_lq3,ext,                               &
                      nx_ext,ny_ext,nz_ext,kdim,                        &
                      var_q,lvl,lvl_coord,units,comment,                &
                      qv_ext,istatus)
  WRITE(6,'(a)')    ' Specific humidity read'
  WRITE(6,'(a,i6)') ' istatus= ',istatus
  WRITE(6,'(a,a)')  ' units= ',units(1)
  IF(istatus /= 1) GO TO 598
!
!-----------------------------------------------------------------------
!
!  Fix-up the q field so that q is the same as the first
!  q above ground.
!
!-----------------------------------------------------------------------
!
  DO j=1,ny_ext
    DO i=1,nx_ext
      qvlast=0.
      DO k=nz_ext,1,-1
        IF(qv_ext(i,j,k) > 1. .OR. qv_ext(i,j,k) < 0.) THEN
          qv_ext(i,j,k)=qvlast
        ELSE
          qvlast=qv_ext(i,j,k)
        END IF
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Integrate temperature from surface to get
!  height of pressure levels.
!  Convert temperature to C for hydrostatic integration
!  Convert specific humidity to mixing ratio
!
!-----------------------------------------------------------------------
!
  DO j=1,ny_ext
    DO i=1,nx_ext
      t_sfc=sfct(i,j)-273.15
      w_sfc=sfcmr(i,j)
      DO k=1,nz_ext
        psnd(k)=p_ext(i,j,k)
        tsnd(k)=t_ext(i,j,k)-273.15
        IF(qv_ext(i,j,k) > 0.) THEN
          mrsnd(k)=1000.*qvtomxr(qv_ext(i,j,k))
        ELSE
          mrsnd(k)=0.
        END IF
      END DO
      CALL getht(nz_ext,nz_ext,                                         &
             tsnd,mrsnd,psnd,                                           &
             terr(i,j),sfcp_pa(i,j),t_sfc,w_sfc,                        &
             htsnd)
      DO k=1,nz_ext
        hgt_ext(i,j,k)=htsnd(k)
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Read u component
!
!-----------------------------------------------------------------------
!
  ext='lw3'
  kdim=nz_ext
  CALL read_laps_data(i4time,dir_lw3,ext,                               &
                      nx_ext,ny_ext,nz_ext,nz_ext,                      &
                      var_u,lvl,lvl_coord,units,comment,                &
                      u_ext,istatus)
  PRINT *, ' U Component read'
  PRINT *, ' istatus= ',istatus
  PRINT *, ' units= ',units(1)
  IF(istatus /= 1) GO TO 598
!
!-----------------------------------------------------------------------
!
!  Read v component
!
!-----------------------------------------------------------------------
!
  ext='lw3'
  kdim=nz_ext
  CALL read_laps_data(i4time,dir_lw3,ext,                               &
                      nx_ext,ny_ext,nz_ext,nz_ext,                      &
                      var_v,lvl,lvl_coord,units,comment,                &
                      v_ext,istatus)
  PRINT *, ' V Component read'
  PRINT *, ' istatus= ',istatus
  PRINT *, ' units= ',units(1)
  IF(istatus /= 1) GO TO 598
!
!-----------------------------------------------------------------------
!
!  Fill qc,qr,qi,qs,qh arrays with missing value.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz_ext
    DO j=1,ny_ext
      DO i=1,nx_ext
        qc_ext(i,j,k)=-999.
        qr_ext(i,j,k)=-999.
        qi_ext(i,j,k)=-999.
        qs_ext(i,j,k)=-999.
        qh_ext(i,j,k)=-999.
      END DO
    END DO
  END DO
!
  istatus=1
  RETURN
!
!-----------------------------------------------------------------------
!
!  Error destination
!
!-----------------------------------------------------------------------
!
  598 CONTINUE
  WRITE(6,'(a)')  ' Error reading last field, returning'
  RETURN
END SUBROUTINE getlaps
!


SUBROUTINE getht(maxlev,nlevel,                                         & 1
           t,w,p_pa,elev,p_sfc_pa,t_sfc,w_sfc,ht)
  IMPLICIT NONE
!
!  Integrate hypsometric equation to get hydrostatic
!  heights from temperatures on pressure levels.
!
!  Input variables:   t  temperature (C)
!                     w  mixing ratio (g/kg)
!                     p_pa pressure in Pascals
!                     elev   elevation of surface in meters
!                     p_sfc_pa  surface pressure in Pascals
!                     t_sfc     surface temperature in C
!                     w_sfc     surface mixing ratio in g/kg
!
!  Output variable:
!                     ht        heights on pressure levels in meters
!
!
!  Keith Brewster Feb, 1994
!
!
!  Input variables
!
  INTEGER :: maxlev,nlevel
  REAL :: t(maxlev),w(maxlev),p_pa(maxlev)
  REAL :: elev,p_sfc_pa,t_sfc,w_sfc
!
!  Output
!
  REAL :: ht(maxlev)
!
!  Parameters
!
  REAL :: g,rd,rdovg
  PARAMETER (g=9.80,rd=287.0)
  PARAMETER (rdovg=(rd/g))
!
!  Functions
!
  REAL :: tctotv
!
!  Misc internal variables
!
  INTEGER :: k,ksfc
  REAL :: tvsfc,tvhi,tvlo,tvbar,thick
!
  DO k=1,nlevel-1
    IF(p_pa(k) < p_sfc_pa) EXIT
  END DO
  ksfc=k
!
  tvsfc=tctotv(t_sfc,w_sfc)
  tvhi =tctotv(  tvbar=0.5*(tvsfc+tvhi)
  ht(ksfc)=elev+rdovg*tvbar*ALOG(p_sfc_pa/p_pa(ksfc))
  IF(ksfc > 1) THEN
    tvlo =tctotv(    tvbar=0.5*(tvsfc+tvlo)
    ht(ksfc-1)=elev+rdovg*tvbar*ALOG(p_sfc_pa/p_pa(ksfc-1))
  END IF
  DO k=ksfc-2,1,-1
    tvlo=tctotv(    tvhi=tctotv(    tvbar=0.5*(tvlo+tvhi)
    ht(k)=ht(k+1)-rdovg*tvbar*ALOG(p_pa(k)/p_pa(k+1))
  END DO
  DO k=ksfc+1,nlevel
    tvlo=tctotv(    tvhi=tctotv(    tvbar=0.5*(tvlo+tvhi)
    ht(k)=ht(k-1)+rdovg*tvbar*ALOG(p_pa(k-1)/p_pa(k))
  END DO
  RETURN
END SUBROUTINE getht
!


  FUNCTION tctotv(tt,ww)
  IMPLICIT NONE
!
!  Virtual Temperature
!
!  Given T in Celcius and mixing ratio in g/kg
!  find the virtual temperature.
!
  REAL :: tctotv,tt,ww
  tctotv=(tt+273.15)*(1.+0.0006*ww)
  RETURN
  END FUNCTION tctotv
!


  FUNCTION qvtomxr(qv)
!
!  This is an approximation to the formula
!
!      qv = w/(1+w)
!
!  Or w = qv + qv*w
!       = qv + qv*(qv + qv*w)
!  So, when q is small, w is approx equal to q
!  so
!     w = qv + qv*(qv + qv*qv)
!
!  Keith Brewster, CAPS
!  Feb 1994
!
!
  IMPLICIT NONE
  REAL :: qvtomxr,qv
!
  qvtomxr=qv + qv*(qv + qv*qv)
!  print *, ' qv, mixr = ',qv,qvtomxr
!
  RETURN
  END FUNCTION qvtomxr