      program create_netcdf

      IMPLICIT NONE

include "/d3/crilmd/LMDZ20091210.trunk/netcdf-4.0.1/include/netcdf.inc"
                 
!--------------------------------------------------------------------------------------
! Lecture du profil initial et des tendances pour le cas FIRE du 14/07/0987 a/c 8:00h
! profil init= z, thetal, qt, u,v,tke                     sur 120 niveaux
! tendances  = z,ugeo,vgeo,wls,dqtdx,dqtdy,dqtdt,thl_rad  sur 120 niveaux
!--------------------------------------------------------------------------------------
! lm=nb niveaux atm., ls= nb niveaux dans sol
! nbt= nb de tendances
! nbstep= nb de tendances
! nbstep2= nb de flux
! nbv= nb variables dans fire.nc
      INTEGER lm,lm2,ls,nbstep,nbstep2,nbt
      PARAMETER (lm=120,lm2=120,ls=1,nbstep=1,nbstep2=1,nbt=8)  ! MPL 04/03/2013
      
      INTEGER k,l,jk,ls2,psurf2,nb_forc,lev_forc
!

!temps
! nbstep = nb de forcages
! nbstep2= nb de flux
      REAL*8 time(nbstep),heure(nbstep)
      REAL*8 time2(nbstep2),hflux(nbstep2)
      REAL*8 time_out(nbstep2)
      INTEGER an,mois,jour
      

!niveaux verticaux
      REAL*8 z_init(lm),z_final(lm)
!     REAL*8 z1(lm),zz(lm2,nbstep),zz2(lm)
      REAL*8 z1(lm),zz2(lm)
      REAL*8 psurf
!     PARAMETER(psurf=98800.)
      REAL*8 pres(lm)

!champs des profils initiaux 
      REAL*8 thetal_init(lm),qt_init(lm),u_init(lm),v_init(lm),tke_init(lm)
      REAL*8 z0,z0t,lati,long,pclay,psand,tg(ls+1),sm(ls+1)
      

!forcages grande echelle
      REAL*8 tend(lm2,nbstep,nbt)

      INTEGER tend_u,tend_v,tend_w,tend_t,tend_q
      INTEGER nudg_u,nudg_v,nudg_w,nudg_t,nudg_q
      

      REAL*8 XLVTT
      PARAMETER(XLVTT=2.5008E+6)

      integer nbvar3d_out
      parameter (nbvar3d_out=27)
      character (len=50), dimension(nbvar3d_out) :: varname3d_out
      character*5 type_sol
      character*4 poub

      integer ierr
      
      integer :: altdimout,timedimout
!,timevarout
      integer :: nout

      integer var3didout(27),toto(lm)

!     Read in data
!     read initial profiles

      open(10,file='profil_init.txt')
      read (10,*) poub
      read (10,*) poub
      do k=1,lm
        read(10,*),z_init(k),thetal_init(k),qt_init(k)&
     ,u_init(k),v_init(k),tke_init(k)
      enddo
      close(10)


!read tendencies
!-------------------------------
! TEND A B C D E
! A=1 on utilise tendances sur U
! B=1 on utilise tendances sur V
! C=1 on utilise tendances sur W
! D=1 on utilise tendances sur theta
! E=1 on utilise tendances sur vapeur eau
!-------------------------------
! NUDG A B C D E
! A= U est nudge avec un temps de relaxation de A
! B= V est nudge avec un temps de relaxation de B
! C= W est nudge avec un temps de relaxation de C
! D= theta est nudge avec un temps de relaxation de D
! E= vapeur d'eau est nudgee avec un temps de relaxation de E
! si A,B,C,D,E=0 => pas de nudging
! si A,B=-1 => nudging avec vent geostrophique fourni
!-------------------------------
      open(11,file='profil_tend.txt')
      read(11,*),poub
      read(11,*),poub
        do k=1,lm2
        read(11,*),tend(k,1,:)
!       print *,'tendances= ',l,k,tend(k,l,:)
        enddo
      close(11)

      do jk=1,lm
!      print *,'jk z pp z_tend p_tend z_final p_final',jk,z(jk),pp(jk),tend2(jk,12,1),tend2(jk,12,2),&
!      tend_final(jk,12,1),tend_final(jk,12,2)
      enddo

!ecriture des resultats dans un fichier nc

         print*,'avant ecriture ok'
 
         call initiate ("fire.nc",z_init,time_out,&
                        altdimout,timedimout,nbstep2,lm)

         print*,'apres initiate'
!        print*,'=============='
         call def_var(nout,"zz","z","m",1,(/altdimout/),var3didout(1),ierr)
         call def_var(nout,"thetal","thetal","K",1,(/altdimout/),var3didout(2),ierr)
         call def_var(nout,"qt","qt","kg/kg",1,(/altdimout/),var3didout(3),ierr)
         call def_var(nout,"u","u","m/s",1,(/altdimout/),var3didout(4),ierr)
         call def_var(nout,"v","v","m/s",1,(/altdimout/),var3didout(5),ierr)
         call def_var(nout,"tke","tke","m2/s2",1,(/altdimout/),var3didout(6),ierr)
         call def_var(nout,"ugeo","ugeo","m/s",2,(/altdimout,timedimout/),var3didout(7),ierr)
         call def_var(nout,"vgeo","vgeo","m/s",2,(/altdimout,timedimout/),var3didout(8),ierr)
         call def_var(nout,"wls","wls","m/s",2,(/altdimout,timedimout/),var3didout(9),ierr)
         call def_var(nout,"dqtdx","dqtdx","m-1",2,(/altdimout,timedimout/),var3didout(10),ierr)
         call def_var(nout,"dqtdy","dqtdy","m-1",2,(/altdimout,timedimout/),var3didout(11),ierr)
         call def_var(nout,"dqtdt","dqtdt","m-1",2,(/altdimout,timedimout/),var3didout(12),ierr)
         call def_var(nout,"thl_rad","thl_rad","K/s",2,(/altdimout,timedimout/),var3didout(13),ierr)
         call def_var(nout,"tend_u","tend_u","-",1,(/altdimout/),var3didout(14),ierr)
         call def_var(nout,"tend_v","tend_v","-",1,(/altdimout/),var3didout(15),ierr)
         call def_var(nout,"tend_w","tend_w","-",1,(/altdimout/),var3didout(16),ierr)
         call def_var(nout,"tend_t","tend_t","-",1,(/altdimout/),var3didout(17),ierr)
         call def_var(nout,"tend_q","tend_q","-",1,(/altdimout/),var3didout(18),ierr)
         call def_var(nout,"nudg_u","nudg_u","-",1,(/altdimout/),var3didout(19),ierr)
         call def_var(nout,"nudg_v","nudg_v","-",1,(/altdimout/),var3didout(20),ierr)
         call def_var(nout,"nudg_w","nudg_w","-",1,(/altdimout/),var3didout(21),ierr)
         call def_var(nout,"nudg_t","nudg_t","-",1,(/altdimout/),var3didout(22),ierr)
         call def_var(nout,"nudg_q","nudg_q","-",1,(/altdimout/),var3didout(23),ierr)

! A FAIRE: ajouter a la suite les scalaires qui sont fournis par profil_init.txt
    
#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(1),z_init)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(1),z_init)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_1"
         endif
 
#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(2),thetal_init)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(2),thetal_init)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_2"
         endif

#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(3),qt_init)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(3),qt_init)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_3"
         endif

#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(4),u_init)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(4),u_init)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_4"
         endif
     
#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(5),v_init)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(5),v_init)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_5"
         endif

#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(6),tke_init)
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(6),tke_init)
#endif     
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_6"
         endif

        do jk=2,nbt
#ifdef NC_DOUBLE
         ierr= NF_PUT_VAR_DOUBLE(nout,var3didout(7+jk-2),tend(:,:,jk))
#else
         ierr= NF_PUT_VAR_REAL(nout,var3didout(7+jk-2),tend(:,:,jk))
#endif
         if(ierr/=NF_NOERR) then
            write(*,*) NF_STRERROR(ierr)
            stop "putvar3d_tendances"
         endif
        enddo

         toto(:)=tend_u
         ierr= NF_PUT_VAR_INT(nout,var3didout(14),toto)
         toto(:)=tend_v
         ierr= NF_PUT_VAR_INT(nout,var3didout(15),toto)
         toto(:)=tend_w
         ierr= NF_PUT_VAR_INT(nout,var3didout(16),toto)
         toto(:)=tend_t
         ierr= NF_PUT_VAR_INT(nout,var3didout(17),toto)
         toto(:)=tend_q
         ierr= NF_PUT_VAR_INT(nout,var3didout(18),toto)
         toto(:)=nudg_u
         ierr= NF_PUT_VAR_INT(nout,var3didout(19),toto)
         toto(:)=nudg_v
         ierr= NF_PUT_VAR_INT(nout,var3didout(20),toto)
         toto(:)=nudg_w
         ierr= NF_PUT_VAR_INT(nout,var3didout(21),toto)
         toto(:)=nudg_t
         ierr= NF_PUT_VAR_INT(nout,var3didout(22),toto)
         toto(:)=nudg_q
         ierr= NF_PUT_VAR_INT(nout,var3didout(23),toto)

ierr=NF_CLOSE(nout)       
       
!
      

       contains
!**********************************************************
Subroutine initiate (filename,z,time,&
                     altdimout,timedimout,&
                     nj,lm)

implicit none

include "/d3/crilmd/LMDZ20091210.trunk/netcdf-4.0.1/include/netcdf.inc"

character (len=*), intent(in):: filename
real*8, dimension(:), intent(in):: z,time
!integer, intent(out):: nout,timevarout
integer, intent(out):: altdimout,timedimout

integer :: altdim,timedim
integer :: newvarid
integer :: nvarid,ierr,i,kloop,nj,lm
!real, dimension(nj), intent(out) :: time

write(*,*) "creating "//trim(adjustl(filename))//'...'
ierr = NF_CREATE(filename,NF_CLOBBER, nout)
if (ierr.NE.NF_NOERR) THEN
   WRITE(*,*)'ERROR: Impossible to create the file.'
   WRITE(*,*) NF_STRERROR(ierr)
   stop ""
endif

!ierr = NF_DEF_DIM (nout, "lat", size(lat), latdimout)
!ierr = NF_DEF_DIM (nout, "lon", size(lon), londimout)
ierr = NF_DEF_DIM (nout, "lev", lm, altdimout)
ierr = NF_DEF_DIM (nout, "time", nj, timedimout)

ierr = NF_ENDDEF(nout)

 call def_var(nout,"lev","lev","m",1,(/altdimout/),nvarid,ierr)
#ifdef NC_DOUBLE
ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,z(1:lm))
#else
ierr = NF_PUT_VAR_REAL (nout,nvarid,z(1:lm))
#endif

 call def_var(nout,"time","time","days since 1987-07-14 08:00:00",1,(/timedimout/),nvarid,ierr)
!ierr = NF_PUT_VAR_INT (nout,timevarout,time(1:nj))
#ifdef NC_DOUBLE
!ierr = NF_PUT_VAR_DOUBLE (nout,timevarout,time(1:nj))
ierr = NF_PUT_VAR_DOUBLE (nout,nvarid,time(1:nj))
#else
!ierr = NF_PUT_VAR_REAL (nout,timevarout,time(1:nj))
ierr = NF_PUT_VAR_REAL (nout,nvarid,time(1:nj))
#endif


end subroutine initiate

!*************************************************************
subroutine def_var(nid,name,title,units,nbdim,dim,nvarid,ierr)

implicit none

include "/d3/crilmd/LMDZ20091210.trunk/netcdf-4.0.1/include/netcdf.inc"

character (len=*) :: title,units,name
integer :: nid,nbdim,nvarid,ierr
integer, dimension(nbdim) :: dim

ierr=NF_REDEF(nid)
#ifdef NC_DOUBLE
ierr = NF_DEF_VAR (nid,adjustl(name),NF_DOUBLE,nbdim,dim,nvarid)
#else
ierr = NF_DEF_VAR (nid,adjustl(name),NF_FLOAT,nbdim,dim,nvarid)
#endif
if(ierr/=NF_NOERR) then
   write(*,*) NF_STRERROR(ierr)
   stop "in def_var"
endif
ierr = NF_PUT_ATT_TEXT (nid, nvarid, "title", len_trim(adjustl(title)),adjustl(title))
if(ierr/=NF_NOERR) then
   write(*,*) NF_STRERROR(ierr)
   stop "in def_var"
endif
ierr = NF_PUT_ATT_TEXT (nid, nvarid, "units", len_trim(adjustl(units)),adjustl(units))
if(ierr/=NF_NOERR) then
   write(*,*) NF_STRERROR(ierr)
   stop "in def_var"
endif
ierr = NF_ENDDEF(nid)

end subroutine def_var

!*************************************************************
subroutine catchaxis(nid,iim,jjm,llm,nj,xc,yc,zc,time,ierr)

!include "/home/distrib/local/netcdf-3.6.1/include/netcdf.inc"
include "/d3/crilmd/LMDZ20091210.trunk/netcdf-4.0.1/include/netcdf.inc"
integer, intent(in) :: nid,iim,jjm,llm,nj
real*8, dimension(iim), intent(out) :: xc
real*8, dimension(jjm), intent(out) :: yc
real*8, dimension(llm), intent(out) :: zc
real*8, dimension(nj), intent(out) :: time
integer, intent(out) :: ierr

integer :: i
integer :: latvar,lonvar,altvar,timevar
integer :: latlen,lonlen,altlen,timelen
integer :: londimin,latdimin,altdimin,timedimin

! Control & lecture on dimensions
! ===============================
   ierr=NF_INQ_DIMID(nid,"yc",latdimin)
   ierr=NF_INQ_VARID(nid,"yc",latvar)
   if (ierr.NE.NF_NOERR) then
      write(*,*) 'ERROR: Field <lat> is missing'
      stop ""  
   endif
   ierr=NF_INQ_DIMLEN(nid,latdimin,latlen)

   ierr=NF_INQ_DIMID(nid,"xc",londimin)
   ierr=NF_INQ_VARID(nid,"xc",lonvar)
   if (ierr.NE.NF_NOERR) then
      write(*,*) 'ERROR: Field <lon> is lacking'
      stop "" 
   endif
  ierr=NF_INQ_DIMLEN(nid,londimin,lonlen)

   ierr=NF_INQ_DIMID(nid,"zc",altdimin)
   ierr=NF_INQ_VARID(nid,"zc",altvar)
   if (ierr.NE.NF_NOERR) then
      write(*,*) 'ERROR: Field <presnivs> is lacking'
      stop ""
   endif
   ierr=NF_INQ_DIMLEN(nid,altdimin,altlen)

!test:on rajoute la dimension temps
   ierr=NF_INQ_DIMID(nid,"time",timedimin)
   ierr=NF_INQ_VARID(nid,"time",timevar)
   if (ierr.NE.NF_NOERR) then
      write(*,*) 'ERROR: Field <time> is lacking'
      stop ""
   endif
   ierr=NF_INQ_DIMLEN(nid,timedimin,timelen)

   if((latlen/=jjm).or.(lonlen/=iim).or.(altlen/=llm)) then
      write(*,*) 'ERROR: Not the good lenght for axis'
      write(*,*) 'longitude: ',lonlen,iim+1
      write(*,*) 'latitude: ',latlen,jjm
      write(*,*) 'presniv: ',altlen,llm
!      stop ""  
   endif

#ifdef NC_DOUBLE
   ierr = NF_GET_VAR_DOUBLE(nid,latvar,yc)
   ierr = NF_GET_VAR_DOUBLE(nid,lonvar,xc)
   ierr = NF_GET_VAR_DOUBLE(nid,altvar,zc)
   ierr = NF_GET_VAR_DOUBLE(nid,timevar,time)
#else
   ierr = NF_GET_VAR_REAL(nid,latvar,yc)
   ierr = NF_GET_VAR_REAL(nid,lonvar,xc)
   ierr = NF_GET_VAR_REAL(nid,altvar,zc)
   ierr = NF_GET_VAR_REAL(nid,timevar,time)
#endif

end subroutine catchaxis


!*************************************************************

!*************************************************************
subroutine giveatt(nid,mean)

Implicit none

!include "/home/distrib/local/netcdf-3.6.1/include/netcdf.inc"
include "/d3/crilmd/LMDZ20091210.trunk/netcdf-4.0.1/include/netcdf.inc"
integer nid,ierr
character (len=100):: att
logical,intent(in) :: mean

ierr = NF_REDEF (nid)
att="GDT 1.3"
ierr= NF_PUT_ATT_TEXT(nid,NF_GLOBAL,'conventions',len_trim(att),att)
if(ierr/=NF_NOERR) then
   write(*,*) NF_STRERROR(ierr)
   stop "in giveatt"
endif

if (mean) then
   att="Mars Climate Database v4.0 - Mean variables"
else
   att="Mars Climate Database v4.0 - Standard deviations"
endif

ierr= NF_PUT_ATT_TEXT(nid,NF_GLOBAL,'history',len_trim(att),att)
if(ierr/=NF_NOERR) then
   write(*,*) NF_STRERROR(ierr)
   stop "in giveatt"
endif

att="LMD-AOPP-ESA-CNES"
ierr= NF_PUT_ATT_TEXT(nid,NF_GLOBAL,'institution',len_trim(att),att)
if(ierr/=NF_NOERR) then
   write(*,*) NF_STRERROR(ierr)
   stop "in giveatt"
endif
ierr = NF_ENDDEF (nid)

end subroutine giveatt

!*************************************************************
subroutine missing_value(nout,nvarid,missing)
IMPLICIT NONE

! useful to watch results with ferret

!include "/home/distrib/local/netcdf-3.6.1/include/netcdf.inc"    
include "/d3/crilmd/LMDZ20091210.trunk/netcdf-4.0.1/include/netcdf.inc"
                                                                                   
INTEGER :: nout,nvarid,ierr
REAL :: missing
ierr = NF_REDEF (nout)
#ifdef NC_DOUBLE
ierr= NF_PUT_ATT_DOUBLE(nout,nvarid,'missing_value',NF_DOUBLE,1,missing)
#else
ierr= NF_PUT_ATT_REAL(nout,nvarid,'missing_value',NF_FLOAT,1,missing)
#endif
iF (ierr.NE.NF_NOERR) THEN
     PRINT*, 'anl_NC: missing value attribution failed'
     WRITE(*,*) 'NF_NOERR', NF_NOERR
     CALL abort
ENDIF                                
ierr=NF_ENDDEF(nout)
return
                                                    
end subroutine missing_value

!****************************************************************
subroutine interp_vertical(nlev_tend,z_tend,p_tend&
  ,nlev_init,z_init,p_init,theta_init,rv_init&
  ,u_init,v_init,temp_init,qv_init&
  ,z_final,p_final,theta_final,rv_final&
  ,u_final,v_final,temp_final,qv_final)
IMPLICIT NONE

!-----------------------------------------------------------------------------------
! Vertical interpolation of initial forcing profiles onto tendencies profiles levels
!-----------------------------------------------------------------------------------
! profil_init.txt nb of levels=nlev_init
! contents:
! z_init, p_init, theta_init, rv_init,u_init,v_init,temp_init,qv_init
!-----------------------------------------------------------------------------------
! profil_tend.txt nb of levels=nlev_tend
!                 nb of forcing times= nbt
! contents:
! z_tend, p_tend, du_tend, dv_tend, dw_tend, dth_tend, drv_tend, dt_tend,dq_tend
!-----------------------------------------------------------------------------------

       INTEGER nlevmax
       PARAMETER (nlevmax=41)
       INTEGER nlev_init,nlev_tend,mxcalc

       REAL*8 z_init(nlev_init), p_init(nlev_init)
       REAL*8 z_tend(nlev_tend), p_tend(nlev_tend)
       REAL*8 theta_init(nlev_init),rv_init(nlev_init)
       REAL*8 u_init(nlev_init),v_init(nlev_init)
       REAL*8 temp_init(nlev_init),qv_init(nlev_init)

       REAL*8 z_final(nlev_tend), p_final(nlev_tend)
       REAL*8 theta_final(nlev_tend),rv_final(nlev_tend)
       REAL*8 u_final(nlev_tend),v_final(nlev_tend)
       REAL*8 Temp_final(nlev_tend),qv_final(nlev_tend)

       INTEGER l,k,k1,k2,kp
       REAL*8 aa,frac,frac1,frac2,fact

!      do l = 1, nlev_tend
!      print *,'p_tend p_init rvi:',l,p_tend(l),p_init(l),rv_init(l)
!      enddo
       do l = 1, nlev_tend
!      print *,'p_tend,p_init',l,p_tend(l),p_init(l)

        if (p_tend(l).ge.p_init(nlev_init)) then

         k1=0
         k2=0

         if (p_tend(l).le.p_init(1)) then

         do k = 1, nlev_init-1
          if (p_tend(l).le.p_init(k) .and. p_tend(l).ge.p_init(k+1)) then
            k1=k
            k2=k+1
!           print *,'p_tend,p_init,l,k,k1,k2',l,k,k1,k2,p_tend(l),p_init(k),p_init(k+1)
!           print *,'++ p_tend p_init rvi:',k1,k2,p_tend(l),p_init(k),rv_init(k)
          endif
         enddo

         if (k1.eq.0 .or. k2.eq.0) then
          write(*,*) 'PB! k1, k2 = ',k1,k2
          write(*,*) 'l,p_tend(l) = ',l,p_tend(l)
         do k = 1, nlev_init-1
          write(*,*) 'k,p_init(k) = ',k,p_init(k)
         enddo
         endif

         frac = (p_init(k2)-p_tend(l))/(p_init(k2)-p_init(k1))
         p_final(l)  = p_tend(l)
         z_final(l)  = z_init(k2)   - frac*(z_init(k2)-z_init(k1))
         theta_final(l) = theta_init(k2)  - frac*(theta_init(k2)-theta_init(k1))
!        print *,'>> l k1 k2 thi1 thi2 thf:',l,k1,k2,theta_init(k1),theta_init(k2),theta_final(l)
         rv_final(l) = rv_init(k2)  - frac*(rv_init(k2)-rv_init(k1))
!        print *,'>> l k1 k2 rvi1 rvi2 rvf:',l,k1,k2,rv_init(k1),rv_init(k2),rv_final(l)
         u_final(l) = u_init(k2)  - frac*(u_init(k2)-u_init(k1))
         v_final(l) = v_init(k2)  - frac*(v_init(k2)-v_init(k1))
         temp_final(l)= temp_init(k2) - frac*(temp_init(k2)-temp_init(k1))
         qv_final(l)= qv_init(k2) - frac*(qv_init(k2)-qv_init(k1))

         else !p_tend>p_init(1)

         k1=1
         k2=2
         frac1 = (p_tend(l)-p_init(k2))/(p_init(k1)-p_init(k2))
         frac2 = (p_tend(l)-p_init(k1))/(p_init(k1)-p_init(k2))
         p_final(l)  = p_tend(l)
         z_final(l)  = frac1*z_init(k1)   - frac2*z_init(k2)
         theta_final(l) = frac1*theta_init(k1)  - frac2*theta_init(k2)
!        print *,'>> l k1 k2 thi1 thi2 thf:',l,k1,k2,theta_init(k1),theta_init(k2),theta_final(l)
         rv_final(l) = frac1*rv_init(k1)  - frac2*rv_init(k2)
!        print *,'<< l k1 k2 rvi1 rvi2 rvf:',l,k1,k2,rv_init(k1),rv_init(k2),rv_final(l)
         u_final(l) = frac1*u_init(k1)  - frac2*u_init(k2)
         v_final(l) = frac1*v_init(k1)  - frac2*v_init(k2)
         temp_final(l)= frac1*temp_init(k1) - frac2*temp_init(k2)
         qv_final(l)= frac1*qv_init(k1) - frac2*qv_init(k2)

         endif ! p_tend<p_init(1)
       endif

       enddo 

!      do l = 1,nlev_init
!      print *,'t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l) ',
!    $        l,t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l)
!      enddo

return
end subroutine interp_vertical

!****************************************************************************

      end 

