c***********************************************************************

	subroutine NLTEdlvr09_ZGRID (n_gcm,  
     @          p_gcm, t_gcm, z_gcm, 
     @          co2vmr_gcm, n2vmr_gcm, covmr_gcm, o3pvmr_gcm ,mmean_gcm,
     @          cpnew_gcm)

c       jul 2011 malv+fgg    First version
c***********************************************************************
                                                
	implicit none                                  
        
	include "dimensions.h"
	include "dimphys.h"
	include 'nltedefs.h'          
	include 'nlte_atm.h'        
	include 'tcr_15um.h'
	include 'nlte_data.h'       
	include 'chimiedata.h'
	include 'conc.h'
                                                
c Arguments 
        integer n_gcm
        real p_gcm(n_gcm), t_gcm(n_gcm)
        real co2vmr_gcm(n_gcm), n2vmr_gcm(n_gcm)
        real covmr_gcm(n_gcm), o3pvmr_gcm(n_gcm)
	real z_gcm(n_gcm)
	real mmean_gcm(n_gcm)
	real cpnew_gcm(n_gcm)

c local variables                               
	integer i, j  , iz                               
!	real  distancia, meanm, gz, Hkm
	real  zmin, zmax, deltazz, deltazzy
	real  nt_gcm(n_gcm)
	real  mmean_nlte(n_gcm),cpnew_nlte(n_gcm)
              
c functions                                     
	external 	hrkday_convert                       
	real 		hrkday_convert                          
                                                
c***********************************************************************


! Define working grid for MZ1D model (NL, ZL, ZMIN) 
! y otro mas fino para M.Curtis (NZ, ZX, ZXMIN = ZMIN

! Para ello hace falta una z de ref del GCM, que voy a suponer la inferior

! Primero, construimos escala z_gcm 

!	z_gcm (1) = zmin_gcm             ! [km]

        !write (*,*) ' iz, p, g, H, z =', 1, p_gcm(1), z_gcm(1)
!	do iz = 2, n_gcm
!	do iz=1,n_gcm
!	   z_gcm(iz)=zlay(iz)/1.e3

!	  meanm = ( co2vmr_gcm(iz)*44. + o3pvmr_gcm(iz)*16. 
!     @               + n2vmr_gcm(iz)*28. + covmr_gcm(iz)*28. ) 
!	  meanm = meanm / n_avog
!	  distancia = ( radio + z_gcm(iz-1) )*1.e5 
!	  gz = gg * masa / ( distancia * distancia ) 
!          Hkm = 0.5*( t_gcm(iz)+t_gcm(iz-1) ) / ( meanm * gz )
!          Hkm = kboltzman * Hkm *1e-5                           ! [km] 
!          z_gcm(iz) = z_gcm(iz-1) - Hkm * log( p_gcm(iz)/p_gcm(iz-1) )

          !write (*,*) iz, p_gcm(iz), gz, Hkm, z_gcm(iz)

!        enddo
! Segundo, definimos los límites del modelo, entre las 2 presiones clave

	! Bottom boundary for NLTE model : Pbottom=2e-2mb=1.974e-5 atm
        jlowerboundary = 1 
	do while ( p_gcm(jlowerboundary) .gt. Pbottom_atm )
	   jlowerboundary = jlowerboundary + 1
	enddo
	zmin = z_gcm(jlowerboundary)
!	write (*,*) ' jlowerboundary, Pmin, zmin =', 
!     @            jlowerboundary, p_gcm(jlowerboundary), zmin

	! Top boundary for NLTE model : Ptop=2e-7mb = 1.974e-5 atm
        jtopboundary = jlowerboundary  
	do while ( p_gcm(jtopboundary) .gt. Ptop_atm ) 
	   jtopboundary = jtopboundary + 1
	enddo
	zmax = z_gcm(jtopboundary)
!	write (*,*) ' jtopboundary, Pmax, zmax =', 
!     @      jtopboundary, p_gcm(jtopboundary),zmax

	deltaz = (zmax-zmin) / (nl-1) 
	do i=1,nl                                      
	    zl(i) = zmin + (i-1) * deltaz
	enddo                                          
!	write (*,*) ' ZL grid:  dz,zmin,zmax ', deltaz, zl(1),zl(nl)
! Creamos el perfil interpolando
	call intersp (    pl,zl,nl,      p_gcm,z_gcm,n_gcm, 2)   ! [atm]
	call intersp (     t,zl,nl,      t_gcm,z_gcm,n_gcm, 1)       
	do i = 1, n_gcm
           nt_gcm(i) = 7.339e+21 * p_gcm(i) / t_gcm(i)    ! [cm-3]	   
	enddo
	call intersp (    nt,zl,nl,     nt_gcm,z_gcm,n_gcm, 2)       
	call intersp (co2vmr,zl,nl, co2vmr_gcm,z_gcm,n_gcm, 1)       
	call intersp ( n2vmr,zl,nl,  n2vmr_gcm,z_gcm,n_gcm, 1)       
	call intersp ( covmr,zl,nl,  covmr_gcm,z_gcm,n_gcm, 1)       
	call intersp (o3pvmr,zl,nl, o3pvmr_gcm,z_gcm,n_gcm, 1)  
	call intersp (mmean_nlte,zl,nl,mmean_gcm,z_gcm,n_gcm,1)
	call intersp (cpnew_nlte,zl,nl,cpnew_gcm,z_gcm,n_gcm,1)
	

	do i = 1, nl

           co2(i) = nt(i) * co2vmr(i)
           n2(i) = nt(i) * n2vmr(i)
           co(i) = nt(i) * covmr(i)
           o3p(i) = nt(i) * o3pvmr(i)

!          	hrkday_factor(i) =  hrkday_convert( t(i),        
!     @        	  co2vmr(i), o3pvmr(i), n2vmr(i), covmr(i) )
	   hrkday_factor(i) = hrkday_convert (mmean_nlte(i),cpnew_nlte(i))

	enddo
                                                
                                              

c  Fine grid for transmittance calculations

	deltazy = (zmax-zmin) / (nzy-1) 
	do i=1,nzy                                      
	    zy(i) = zmin + (i-1) * deltazy	             
	enddo                                          
!	write (*,*) ' ZY grid:  nzy,dzy,zmin,zmax ', 
!     @         nzy, deltazy, zy(1),zy(nzy)

	call intersp (    py,zy,nzy,      p_gcm,z_gcm,n_gcm, 2)   ! [atm]
	call intersp (    ty,zy,nzy,      t_gcm,z_gcm,n_gcm, 1)       
	call intersp (   nty,zy,nzy,     nt_gcm,z_gcm,n_gcm, 2)       

	call intersp (  co2y,zy,nzy,   co2vmr_gcm,z_gcm,n_gcm, 1) 
	do i=1,nzy 
	   co2y(i) = co2y(i) * nty(i)
	enddo


  

c end                                           
	return                                         
	end