module modairsco
  use kinds
  implicit none
  integer(i_kind) :: levco(10)=(/1, 20, 45, 56, 63, 70, 81, 89, 93, 100/)
  integer(i_kind), parameter :: nlevco=10
  integer(i_kind) :: lev_1(nlevco)
  integer(i_kind), parameter :: nlairs=100,ncor=9,ncorp=10
  integer nllpointmin
!  integer(i_kind) :: nlev_1=10
! this converted from integral of ppv to molecules/ccm2
! we want to have ppmv for our uses in gsi times dp(mb)
! want to have a number on layers of order 100
  !real(r_kind),parameter :: fac=1.0e4*6.023e23/28.97/9.8/1000.
! make it with ppmv input
!  real(r_kind),parameter :: fac=fac*1.e-6
! converts ppmv to mol/ccm2
! for intco and stpco want to scale from ppmv to mol/ccm2*1.e-15
  logical :: logicalairsco=.false.
  logical :: ltropomico=.false.
!  logical :: debugco=.true.
  logical :: debugco=.false.
  logical :: usewdpij=.false.
!  real(r_kind),parameter :: fac=1.0e4*6.023e23/28.97/9.8/1000.*1.e-6
  real(r_kind) :: aknorm=-1.
  real(r_kind),parameter :: fac=1.0e4*6.023e23/28.97/9.8/1000. ! put back to ppv
  real(r_single),dimension(nlairs,12) :: guess_profile_nh,guess_profile_sh
  real(r_single) :: totprt(nlairs),dpairs(nlairs)
  real(r_single), dimension(100) :: pobsairsco=(/0.0161, 0.0384, 0.0769, 0.1370, 0.2244, 0.3454, 0.5064, 0.7140, &
      0.9753,&
      1.2972,  1.6872,  2.1526,  2.7009,  3.3398,  4.0770,  4.9204,  5.8776, 6.9567,  8.1655,&
      9.5119, 11.0038, 12.6492, 14.4559, 16.4318, 18.5847, 20.9224, 23.4526, 26.1829, 29.1210,&
      32.2744, 35.6505, 39.2566, 43.1001, 47.1882, 51.5278, 56.1260, 60.9895, 66.1253, 71.5398,&
      77.2396, 83.2310, 89.5204, 96.1138,103.0172,110.2366,117.7775,125.6456,133.8462,142.3848, &
      151.2664,160.4959,170.0784,180.0183,190.3203,200.9887,212.0277,223.4415,235.2338,247.4085,&
      259.9691,272.9191,286.2617,300.0000,314.1369,328.6753,343.6176,358.9665,374.7241,390.8926,&
      407.4738,424.4698,441.8819,459.7118,477.9607,496.6298,515.7200,535.2322,555.1669,575.5248,&
      596.3062,617.5112,639.1398,661.1920,683.6673,706.5654,729.8857,753.6275,777.7897,802.3714,&
      827.3713,852.7880,878.6201,904.8659,931.5236,958.5911,986.0666,1013.9476,1042.2319,1070.917,1100.0000/)
  real(r_single), dimension(101) :: pobs2=(/0.005,0.0161, 0.0384, 0.0769, 0.1370, 0.2244, 0.3454, 0.5064, 0.7140, &
      0.9753,&
      1.2972,  1.6872,  2.1526,  2.7009,  3.3398,  4.0770,  4.9204,  5.8776, 6.9567,  8.1655,&
      9.5119, 11.0038, 12.6492, 14.4559, 16.4318, 18.5847, 20.9224, 23.4526, 26.1829, 29.1210,&
      32.2744, 35.6505, 39.2566, 43.1001, 47.1882, 51.5278, 56.1260, 60.9895, 66.1253, 71.5398,&
      77.2396, 83.2310, 89.5204, 96.1138,103.0172,110.2366,117.7775,125.6456,133.8462,142.3848, &
      151.2664,160.4959,170.0784,180.0183,190.3203,200.9887,212.0277,223.4415,235.2338,247.4085,&
      259.9691,272.9191,286.2617,300.0000,314.1369,328.6753,343.6176,358.9665,374.7241,390.8926,&
      407.4738,424.4698,441.8819,459.7118,477.9607,496.6298,515.7200,535.2322,555.1669,575.5248,&
      596.3062,617.5112,639.1398,661.1920,683.6673,706.5654,729.8857,753.6275,777.7897,802.3714,&
      827.3713,852.7880,878.6201,904.8659,931.5236,958.5911,986.0666,1013.9476,1042.2319,1070.917,1100.0000/)
  real(r_single),parameter :: mw_d= 28.9644              ! gm/mole dry air
  real(r_single),parameter :: mw_w=18.0151              ! gm/mole water
  real(r_single),parameter :: g_std=980.665              ! PT99, standard gravity, cm/s^2
  real(r_single),parameter :: eps=mw_w/mw_d            ! ratio of wet to dry air
  real(r_single),parameter :: Navog=6.0221367E+23        ! avogadro's number
  real(r_single),parameter :: qx= 1000.0*Navog/(mw_d*g_std)
  contains
  subroutine getakmatrix(psrf,nlev,nlev_1,lev_1,ak,ncord,&
     tranakfinv)
  use kinds
  implicit none
  real(r_single) :: blmult
  
  real(r_single) :: psrf,coapcd(nlairs),tots(ncord),cotslb(ncord),h2oslb(ncord),coapslb(ncord)
  integer(i_kind) :: ip,il,l1,l2,style,htp,nlev_1,lev_1(nlev_1),nlev,ncord
  logical :: hbt
  real(r_double) :: finv(ncord,nlev),finex(ncord,nlev)
  real(r_double) :: akfinv(ncord,nlev),tranakfinv(nlev,nlev),ak(ncord,ncord)
  real(r_kind) :: coapcda(nlev)
  real(r_kind) :: comrt(ncord),coapmr(ncord)
  style=1
  htp=1
  hbt=.true.
! nlev_1 is variable from changlev not a constant as is nlev and lev_1
! ncord should be nlev_1-1 here
!  write(6,*)'call calc'
!  call flush(6)
  call calc_finv(style,nlev_1,lev_1,nlev,htp,hbt,pobsairsco,nlairs,finv,finex,ncord)
!  write(6,*)'did call calc'
!  call flush(6)
    ! ===================================================
    ! the following is the more correct way of convolving 
    !   
    ! convert the profile to coarse slab values using the 
    ! effective averaging kernel Equations 13, 16, and 17  
    ! of MB08


    ! ==============================================
    ! convolve the "truth" with the averaging kernel 
    ! the equation is : 
    ! ln( co' ) = ln(co_apriori) + F*A*Finv * ln(cotrue / co_apriori) 
    !   
    ! Note: # Computes array elements by multiplying the columns of the 
    ! first array by the rows of the second array. The second array must 
    ! have the same number of columns as the first array has rows. The 
    ! resulting array has the same number of columns as the first array 
    ! and the same number of rows as the second array.    
!  matrix by itself form
!  ak now ndimco x ndimco, finv is now ncord,nlev
   akfinv=matmul(ak,finv)
!   matrix by itself form
   tranakfinv=matmul(transpose(finex),akfinv)
  return
  end subroutine getakmatrix
!PRO CALC_FINV, style, nlev_1, lev_1, numlev, htp, hbt, Pobs, Finv, $
!               fine=finex
! 
! NAME: calc_finv
! PURPOSE: calculate the inverse of the trapezoid matrix 
! 
! INPUT: 
!
!  Name          Description 
! ------    ----------------------
! style     1 for Moore-Penrose pseudoinverse
!           2 for slab averages
! nlev_1    number of coarse levels (topography corrected)
! lev_1     coarse level boundaries 
! numlev    number of levels (RTA grid)  
! htp       halftop
! hbt       halfbot
! Pobs      Pressure of RTA grid

! OUTPUT  

! Finv      Inverse matrix of the trapezoids
! 
! nlev_1, lev_1, numlev all come from changlev and can vary
! ncor is now nlev_1-1
  subroutine calc_finv(style,nlev_1, lev_1, numlev, htp, hbt, Pobs, nlairs,Finv,finex,ncor)
  use kinds
  implicit none
! input
  integer(i_kind),intent(in) :: nlev_1,numlev,style,htp,nlairs,ncor
  logical, intent(in) :: hbt
  real(r_single),intent(in) :: pobs(nlairs)
  integer(i_kind) :: lev_1(nlev_1)
! output
  real(r_double),intent(out) :: finex(ncor,numlev)
!  real(r_double),intent(out) :: finex(nlev_1-1,numlev)
  real(r_double),intent(out) :: finv(ncor,numlev)
  real(r_double) :: finvnew(ncor,numlev)
! local
  real(r_double) :: finea(ncor,numlev),fftr(ncor,ncor),finv1(ncor,ncor)
  real(r_double) :: total
!  real(r_single) :: SLBVAL(nlev_1-1),fine(numlev)
  real(r_single) :: SLBVAL(ncor),fine(numlev)
  real(r_single) :: p1100=1100.
  integer(i_kind) ::status,n,iavg,ii,il,ndim,i,j
  ndim = nlev_1-1
!  write(6,*)'style',style
!  call flush(6)
  select case(style)
    case (1)
      ! calculate trapezoids the trapezoids
!     finea = FLTARR(ndim, numlev)
      do il=1,ndim
!          slbval = FLTARR(ndim)
        slbval=0.0
        slbval(il) = 1.0
!          fine = FLTARR(numlev)
        call slb2fin(nlev_1, lev_1, slbval, htp, hbt,p1100, Pobs, nlairs,fine,numlev )        
        finea(il,:) = fine    !/TOTAL(fine)
      end do

      finex = finea

      
      ! calculate the interpolation matrix 
      fftr=matmul(finea,transpose(finea))
   
      
      status = 1
      finv1=0.0
      call la_invert(fftr,finv1,ncor,status)

    
     ! =============================================================== 
     ! this is the inverse of the interpolation functions (trapezoids)
     ! retrieval basis
      finv=matmul(finv1,finea)
  
!     END
    CASE (2)
!   2L: BEGIN
!     calculate trapezoids the trapezoids
!     finea = FLTARR(ndim, numlev)
      DO IL=1,NDIM
!       slbval = FLTARR(ndim)
        slbval(il) = 1.0
!       fine = FLTARR(numlev)
        call slb2fin(nlev_1, lev_1, slbval, htp, hbt,p1100, Pobs, nlairs,fine,numlev)
        finea(il,:) = fine
      end do
      finex = finea
!     finv = DBLARR(ndim,numlev)
      DO IL=1,NDIM
        N = lev_1(il+1) - lev_1(il) + 1
!       iavg = lev_1(il) + LINDGEN(N) - 1L ! offset for IDL indx
        finv(il,iavg) = 1.0/float(N)  ! all elements are normalized
       
!       ENDFOR
      end do
      do il=1,ndim
        total=0
        do ii=1,numlev
          total=total+finea(il,ii)
        enddo
!       finea(il,:) = finea[il,:]/TOTAL(finea[il,:])
        finea(il,:) = finea(il,:)/TOTAL
      end do
              
      finv = finea  ! just take the transpose
     
!     finv = SMOOTH(finv,2)
!     stop
    CASE DEFAULT
      write(6,*)'default case'
      call flush(6)
      STOP
    END SELECt
   return
END  subroutine calc_finv
subroutine changlev(Psurf, Pstd, nlairs,numlev, nlev_0, lev_0, nlev_1, lev_1)
  use kinds
  implicit none
  integer(i_kind),intent(in) :: nlev_0,lev_0(nlev_0),nlairs
  integer(i_kind) :: nlev,numlev
  integer(i_kind),intent(out) :: nlev_1
  integer(i_kind) lev(nlev_0),L,numlevnew,numlevold
  integer(i_kind),intent(out) :: lev_1(nlev_0)
  real(r_single),intent(in) :: pstd(nlairs),psurf
  real(r_single) :: ptop,pbot,thinlev

  thinlev = 50.0

  numlevold = numlev
  numlevnew = numlev

  ptop=100.
  pbot=1100.
  call lsurface( numlevnew, Pstd, nlairs,Psurf, ptop, pbot)
  
  numlev = numlevnew
  nlev = nlev_0
  lev  = lev_0 
  do L=1,nlev_0
    if ( lev(L) > numlevnew .and.l > 1) then 
      if (lev(l-1) == numlevnew) then 
        nlev = nlev - 1
      else 
        if ( Pstd(numlevnew) - Pstd(lev(l-1)) < thinlev ) then 
          nlev = nlev - 1
          lev(L-1) = numlevnew
        endif
      endif
      lev(L) = numlevnew
    endif
  end do
  nlev_1=nlev
  lev_1  = lev(1:nlev)
  return
end subroutine changlev
! input:
!   numlev = # of levels (e.g., 100)
!   pres = Pobs(0:numlev-1)
!   Psurf = surface pressure
!   Plow = lower limit of pressure for numlev
!
! output:
!   numlev = bottom level index in Pobs(1) system.
!     NOTE: for IDL numlev-1 is the bottom level


!pro lsurface, numlev, pres, Psurf, Plow, Phigh
subroutine lsurface(numlev, pres, nlairs,Psurf, Plow, Phigh)
  use kinds
  implicit none
  integer(i_kind),intent(in) :: nlairs
  real(r_single),intent(in) :: pres(nlairs),psurf,plow
  real(r_single),intent(inout) :: phigh
  integer(i_kind),intent(out) :: numlev
  integer(i_kind) :: L,llsurface
  phigh = 1100.
  if ( psurf > plow .and. psurf <= phigh ) then 
    do L = numlev, 1, -1 
      if ( psurf >= pres(L-1)+5.0 ) then 
        llsurface = L
! will be 1 bigger than idl version
        go to 100
      endif
    end do
  endif
  if(psurf <= plow) then 
!    llsurface = 0
    llsurface = 1
    do L = 1, numlev
      if(psurf < pres(L)) llsurface = L
    end do
  else 
!    llsurface = numlev-1
    llsurface = numlev
    do L = numlev,1,-1 
      if(psurf > pres(L)+5.0) llsurface = L
    end do
  endif

100 continue
  numlev=llsurface
!  numlev=  llsurface + 1
!  numlev=numlev -1 ! make idl index again

end subroutine lsurface
subroutine la_invert(ain,aout,mn,status)
! assume square matrixs ain(mn,mn),aout(mn,mn)
  use kinds
  implicit none
  real(r_double) :: ain(mn,mn),aout(mn,mn)
  integer(i_kind) :: status,mn
  integer(i_kind) :: ipiv(mn)
  integer(i_kind) :: info,lwork,i,j,jj
  real(r_double) :: aprod(mn,mn),work(mn+100000),work2(mn+100000)
  aout=ain
  call dgetrf(mn,mn,aout,mn,ipiv,info)
  info=0
  lwork=100000+mn
  work=0.0
  lwork=-1
  call dgetri(mn,aout,mn,ipiv,work2,lwork,info)
  lwork=100000+mn
  call dgetri(mn,aout,mn,ipiv,work2,lwork,info)
  return
end subroutine la_invert
! same function as the FORTRAN version.  Builds trapezoids.
! numslab = number of levels in Lslab
! Lslab   = index values for trapezoid hinge points, starting at 1
! tmpslab = amplitude of trapezoids
! usehalftop = 1 for trapezoid, 0=wedge
! Psurf = pressure at surface
! Pres  = pressure values at index points
! tmpfine = amplitude on fine levels
! tmpfile(L) = tmpslab_i*F_i(L)


subroutine slb2fin( numslab, Lslab, tmpslab, usehalftop, usehalfbot,&
    Psurf, Pres, nlairs,tmpfine,nfine)
    use kinds
    implicit none   
!     tslab = fltarr(numslab)
!   input variables
    integer(i_kind),intent(in) :: numslab,nlairs,nfine
    integer(i_kind),intent(in) ::lslab(numslab),usehalftop
    real(r_single),intent(in) :: pres(nlairs),psurf,tmpslab(numslab)
    logical,intent(in) :: usehalfbot
!   output variable
    real(r_single), intent(out) :: tmpfine(nfine)
!   local variables
    real(r_single) :: tslab(numslab),ztmpdn
    real(r_single) :: ztmpup,zprlndn,zsloe,zprlnup,zslope
    integer(i_kind) :: ldn,lup,L,n
  

!   ------------------------------------------------------
!   Interpolate temperatures ( linear in log of pressure )
!   ------------------------------------------------------

    if(usehalftop > 0) then
      tslab(1) = 0.5 * tmpslab(1)
    else
      tslab(1) = tmpslab(1)
    endif

    do n=2,numslab-1
      tslab(n) = 0.5 * ( tmpslab(n) + tmpslab(n-1) )
    end do

    if ( usehalfbot ) then
      tslab(numslab) = 0.5 * tmpslab(numslab-1)
    else
      tslab(numslab) = tmpslab(numslab-1)
    endif
    ldn          = lslab(1)
    ztmpdn       = tslab(1)
    zprlndn      = log (pres(ldn) )
    do n=1,numslab-1
      lup        = ldn
      ztmpup     = ztmpdn
      zprlnup    = zprlndn
!      ldn        = lslab(n+1) - 1
      ldn        = lslab(n+1) 
      ztmpdn     = tslab(n+1)
!     zprlndn    = alog(Pres(ldn))
      zprlndn    = log(Pres(ldn))
      zslope     = ( ztmpdn - ztmpup ) / ( zprlndn - zprlnup )
        
      do L=lup,ldn-1
        tmpfine(L) = ztmpup + zslope * ( log(Pres(L)) - zprlnup )
      end do
      
     lup = L
  
    end do
      
    tmpfine(ldn) = ztmpdn
    return
end subroutine slb2fin

    subroutine readcoguessprofile
    implicit none
    integer(i_kind) :: l,i
    real(r_kind) :: pressure(nlairs)
    open(30,file='V6_CO_Initial_Guess_Profiles_NH.csv',form='formatted')
    do l=1,nlairs
      read(30,*)pressure(l),(guess_profile_nh(l,i),i=1,12)
    end do
    open(30,file='V6_CO_Initial_Guess_Profiles_SH.csv',form='formatted')
    do l=1,nlairs
      read(30,*)pressure(l),(guess_profile_sh(l,i),i=1,12)
    end do
    do l=1,nlairs
      totprt(l)=qx*(pobs2(l+1)-pobs2(l))
      dpairs(l)=pobs2(l+1)-pobs2(l)
    end do
    nllpointmin=-999
    
    end subroutine readcoguessprofile


end module modairsco
