module mwoa

use mkinds
use mdates
use mncio
use mutil
use mseawater
use mgdem
use mgrid
use mhycom
use miotsis
use mpobs

implicit none

contains
subroutine woa2tsis_profiles(cdtg,cfname)
implicit none
character(22)             :: cfname
character(10)             :: cdtg,cdtgl
integer                   :: ncid,nprof,nlvl,vid,fid
real(r4),allocatable      :: lon(:),lat(:),depth(:)
real(r4),allocatable      :: ct(:,:,:),cs(:,:,:)
integer(i4),parameter     :: ntzlevs=102,cnx=360,cny=180
real(r4)                  :: tzlevels(ntzlevs)
integer                   :: i,j,ncount,k,kk,stat,ngdem,fcount,ii,jj,kmx,l,nlp
character(:),allocatable  :: CFileName
character(:), allocatable :: vname
integer(i4)               :: vtype
integer(i4)               :: vdims3(3),vdims2(2),vdim(1)
real(r4)                  :: date_offset,wday,wsday
logical(bl)               :: lexist
integer,parameter         :: nprofmx=10000
real(r4),allocatable      :: plon(:),plat(:),ptemp(:,:),psaln(:,:),ptime(:),pdens(:,:)
integer(i4),allocatable   :: indx(:,:),akmx(:,:),prof_kmx(:),prof_valid(:),prof_num(:)
real(r4),allocatable      :: pl(:),rl(:),sl(:),tl(:),ul(:),vl(:),dp(:),terr(:),serr(:),derr(:),thkerr(:),tzerr(:),szerr(:)
real(r4),allocatable      :: plt(:),rlt(:),slt(:),tlt(:),iplt(:),tzt1(:),tzt2(:)
character(:),allocatable  :: var
real(r4),parameter        :: spval=2.00**100
real(r4)                  :: rho_min,tsig,q
! error dp or d
real(r4),parameter        :: sigdmin  = 0.005
real(r4),parameter        :: sigdpmin = 0.5
real(r4),parameter        :: factdp   = 0.02
real(r4),parameter        :: fact_tserr=0.5
real(r4)                  :: fact_tserr_tau
real(r4),parameter        :: argo_temp_err=0.01, argo_sal_err=0.005
real(r4)                  :: il,jl

!model info
character(len=20)         :: model,fcst_file,anl_file
logical(bl)               :: lglb,fcst_file_out
character(len=5)          :: pfx
NAMELIST/model_info/model,sigver,lglb,fcst_file,anl_file,fcst_file_out,pfx

! reject counts
integer(i4)               :: ndate_time_qc_fail,ninv_qc_fail,nclim_qc_fail
real(r4)                  :: ipr,jpr,min_rho

date_offset=ndays(01,01,1950,12,31,1900)

! set target z levels


call getarg(1,cdtg)
call getarg(2,cfname)
!call CDATE2WNDAY(WDAY,cdtg)
call read_profile_params()
call initialize_analysis_grid()
call read_obs_location()
call read_obs_window()
call initialize_pobs_vector()

OPEN(21,FILE='tsis.nlist')
READ(21,NML=model_info)
CLOSE(21)

call nciopn(cfname,fid)
call ncioin(cfname,fid,"nprof",ncount)
call nciorv(cfname,fid,"depth",tzlevels)

if(ncount>1) then

allocate(indx(nx,ny))
allocate(ptime(ncount))
allocate(plon(ncount))
allocate(plat(ncount))
allocate(ptemp(ncount,ntzlevs))
allocate(psaln(ncount,ntzlevs))
allocate(pdens(ncount,ntzlevs))
allocate(tzerr(ntzlevs))
allocate(szerr(ntzlevs))
allocate(tzt1(ntzlevs))
allocate(tzt2(ntzlevs))
allocate(pl(nz+1),rl(nz),sl(nz),tl(nz),ul(nz),vl(nz),dp(nz),derr(nz),thkerr(nz),terr(nz),serr(nz))
allocate(plt(nz+1),rlt(nz),slt(nz),tlt(nz),iplt(nz))

call nciorv(cfname,fid,"plon",plon)
call nciorv(cfname,fid,"plat",plat)
call nciorv(cfname,fid,"temp",ptemp)
call nciorv(cfname,fid,"terr",tzerr)
call nciorv(cfname,fid,"saln",psaln)
call nciorv(cfname,fid,"serr",szerr)
call nciocl(cfname,fid)



! get std deviations

fcount=ncount
   ncount=0
   do i=1,fcount
           if(lglb) then
             call lon_lat_idx(mlon,mlat,nx,ny,plon(i),plat(i),lglb,il,jl)
             ii=int(il)
             jj=int(jl)
           else
             ii=int(ridx(mlon(:,1),plon(i)))
             jj=int(ridx(mlat(1,:),plat(i)))
           endif
           if(mdepth(ii,jj)<spval) then
           indx(ii,jj)=indx(ii,jj)+1
           if(indx(ii,jj)==1 ) then
              kmx=ntzlevs 
           do k=1,ntzlevs
            if(psaln(i,k)==-999.0) then
             kmx=k-1
             exit
            endif
           enddo

            !write(*,*) "Profile No: ", ncount,plon(i),plat(i)
            !write(*,*) "****************************************************************"
            do k=1,ntzlevs
             if(k<=kmx) then
              pdens(i,k)= sig(ptemp(i,k),psaln(i,k),sigver)
             else
              pdens(i,k)=-999.0       
             endif

            !write(*,'(I3.3,4F16.4)') k,tzlevels(k),ptemp(i,k),psaln(i,k),pdens(i,k)
            enddo
           ! layerize the profile
            pl=0.0
            thkmin=1.0
            mlayer=isosigma(ii,jj,:) ! needed for vsigma
            !print *, i, mdepth(ii,jj),tzlevels(kmx)
           call layerize(kmx,tzlevels(1:kmx),ptemp(i,1:kmx),psaln(i,1:kmx),pdens(i,1:kmx),tzerr,szerr,min(tzlevels(kmx),mdepth(ii,jj)),&
           nz,pl,rl,sl,tl,terr,serr,thbase,thkmin,nhybrd,nsigma,dp00,dp00x,dp00f,ds00,ds00x,ds00f,dp0k,ds0k(1:nsigma),ilspec,mlayer,thkmin,sigver)
           pl(1)=0.0

           ! find the last layer with info
            do k=1,nz-1
            if(abs(tl(K+1)-tl(k))<bot_check_tol .and. abs(sl(K+1)-sl(k))<bot_check_tol .and. abs(pl(K+1)-pl(k))<bot_check_tol )then
                kmx=k
            exit
            endif
           enddo

!           diag info
           

           do k=1,nz
               if(k<=kmx) then
                 derr(k)=min(abs(dsig(tl(k),sl(k),terr(k),serr(k),sigver)),max(abs(mlayer(k)-rl(k)),sigdmin))
                 derr(k)=max(sigdmin,derr(k))
                 thkerr(k)=max(sigdpmin, (factdp*pl(k) -(factdp*pl(k)-sigdpmin)*(derr(k)-sigdmin)/derr(k)))
                 thkerr(k)=min(100.0,thkerr(k))
               else
                 derr(k)=sigdmin
                 thkerr(k)=sigdpmin

                 !write(*,'(I5,8F16.4)') k,pl(k),tl(k),terr(k),sl(k),serr(k),rl(k),derr(k),thkerr(k)

                endif
             enddo

           

           ncount=ncount+1
           !write(*,*) "Profile No: ", ncount
           !write(*,*) "****************************************************************"
           do k=1,nz
           !if(k<kmx) write(*,'(I5,8F16.4)') k,pl(k),tl(k),terr(k),sl(k),serr(k),rl(k),derr(k),thkerr(k)
             work_pobs(1,k,ncount)%lon=plon(i)
             work_pobs(1,k,ncount)%lat=plat(i)
             work_pobs(1,k,ncount)%val=tl(k)
             work_pobs(1,k,ncount)%err=terr(k)
             work_pobs(1,k,ncount)%gridi=ii
             work_pobs(1,k,ncount)%gridj=jj
             if(k<kmx) then
                work_pobs(1,k,ncount)%id=1
             else if(k>=kmx) then
                work_pobs(1,k,ncount)%id=0
             endif

             work_pobs(2,k,ncount)%lon=plon(i)
             work_pobs(2,k,ncount)%lat=plat(i)
             work_pobs(2,k,ncount)%val=sl(k)
             work_pobs(2,k,ncount)%err=serr(k)
             work_pobs(2,k,ncount)%gridi=ii
             work_pobs(2,k,ncount)%gridj=jj
             if(k<kmx) then
                work_pobs(2,k,ncount)%id=1
             else if(k>=kmx) then
                work_pobs(2,k,ncount)%id=0
             endif

             work_pobs(3,k,ncount)%lon=plon(i)
             work_pobs(3,k,ncount)%lat=plat(i)
             work_pobs(3,k,ncount)%val=rl(k)
             work_pobs(3,k,ncount)%err=derr(k)
             work_pobs(3,k,ncount)%gridi=ii
             work_pobs(3,k,ncount)%gridj=jj
             if(k<kmx) then
                work_pobs(3,k,ncount)%id=1
             else if(k>=kmx) then
                work_pobs(3,k,ncount)%id=0
             endif

             work_pobs(4,k,ncount)%lon=plon(i)
             work_pobs(4,k,ncount)%lat=plat(i)
             work_pobs(4,k,ncount)%val=pl(k)!*9806.
             work_pobs(4,k,ncount)%err=thkerr(k)!*9806.
             work_pobs(4,k,ncount)%gridi=ii
             work_pobs(4,k,ncount)%gridj=jj
             if(k<kmx) then
                work_pobs(4,k,ncount)%id=1
             else if(k>=kmx) then
                work_pobs(4,k,ncount)%id=0
             endif

           enddo  !k
            endif ! land/sea
           endif ! uniq
!      endif   ! prof_valid
    enddo     ! nprof
   
!   enddo
! endif ! super

! tsis_profile_model_grid
allocate(profile_obs(nobs_typ,nz,ncount))
profile_obs=work_pobs(:,:,1:ncount)
endif
call write_tsis_profiles(cdtg,ncount,nz,nobs_typ,mlon(:,1),mlat(1,:),nx,ny,pfx,1)
write(*,*) "Total no of profiles: ", ncount

!print *, fcount,ncount
end subroutine


end module