program ERAi
use netcdf
implicit none

character(len=132)  ::  modfile,modfileout,modfile2,modfileout2,modfile3
integer  ::  ix,iy,iz,it,id,ix2,iy2
integer,parameter  ::  tempmax=38,lonera=480,latera=241,altmax2=37,daymax=12,altmax=40
integer,parameter  ::  lonmax=160,latmax=80
real,parameter  ::  undef=-9999.
real,parameter :: Aw=-6096.9385, Bw=16.635794,Cw=-2.711193e-02,Dw=1.673952e-05,Ew=2.433502
real,parameter :: Ai=-6024.5282,Bi=24.7219,Ci=1.0613868e-02,Di=-1.3198825e-05,Ei=-0.49382577
real :: A,B,C,D,E,temptmp
real  ::  min
real,dimension(tempmax+1)  ::  tempmod
real,dimension(altmax2)  ::  presnivs
!real,dimension(lonera,latera,altmax,daymax)  ::  temp,rhum


real,dimension(altmax2,tempmax,daymax)  ::  vap,indvap
real,dimension(lonera,latera,altmax2,daymax)  ::  temp,rhum


real,dimension(altmax)  ::  pres_goccp


modfile='/bdd/ERAI/NETCDF/GLOBAL_075/1xmonthly/AN_PL/2010/r.2010.apmei.GLOBAL_075.nc'
modfile2='/bdd/ERAI/NETCDF/GLOBAL_075/1xmonthly/AN_PL/2010/ta.2010.apmei.GLOBAL_075.nc'

modfileout='/homedata/gcesana/DEPOL_restore/DEPOL/ERAI_vap.nc'


tempmod=(/-120.,-90.,-87.,-84.,-81.,-78.,-75.,-72.,-69.,-66.,-63.,-60.,-57.,-54.,  &
          -51.,-48.,-45.,-42.,-39.,-36.,-33.,-30.,-27.,-24.,-21.,-18.,-15.,-12.,   &
          -9.,-6.,-3.,0.,3.,6.,9.,12.,15.,18.,60. /)

tempmod=tempmod+273.15

presnivs = (/1000., 975., 950., 925., 900., 875., 850., 825., 800., 775., 750.,    &
            700., 650., 600., 550., 500., 450., 400., 350., 300., 250., 225., 200., &
            175., 150., 125., 100., 70., 50., 30., 20., 10., 7., 5., 3., 2., 1./)


!allocate(rhum(lonera,latera,altmax2,daymax),temp(lonera,latera,altmax2,daymax))
!allocate(vap(altmax2,tempmax,daymax),indvap(altmax2,tempmax,daymax))



print *, 'lecture fichiers modele'
call rdnc4(trim(modfile2),lonera,latera,altmax2,daymax,temp,'ta')
call rdnc4(trim(modfile),lonera,latera,altmax2,daymax,rhum,'r')


do id=1,daymax
print *, 'time = ',id
do iz=1,altmax2
do iy=1,latera
do ix=1,lonera


A=Ai; B=Bi; C=Ci; D=Di; E=Ei
    do it=1,17 !! ice
    if( (temp(ix,iy,iz,id).ge.tempmod(it)).and.  &
        (temp(ix,iy,iz,id).lt.tempmod(it+1)).and.(temp(ix,iy,iz,id).ne.undef) &
       .and.(rhum(ix,iy,iz,id).gt.0.1) )then


temptmp=temp(ix,iy,iz,id)
vap(it,iz,id)=vap(it,iz,id)+rhum(ix,iy,iz,id)*  &
exp(A*(1/temptmp)+B+C*temptmp+D*temptmp**2+E*log(temptmp))/100
indvap(it,iz,id)=indvap(it,iz,id)+1.

    endif
    enddo



A=Ai; B=Bi; C=Ci; D=Di; E=Ei
    do it=18,31  !! mixte
    if( (temp(ix,iy,iz,id).ge.tempmod(it)).and.  &
        (temp(ix,iy,iz,id).lt.tempmod(it+1)).and.(temp(ix,iy,iz,id).ne.undef) &
       .and.(rhum(ix,iy,iz,id).gt.0.1) )then

temptmp=temp(ix,iy,iz,id)
vap(it,iz,id)=vap(it,iz,id)+rhum(ix,iy,iz,id)*  &
exp(A*(1/temptmp)+B+C*temptmp+D*temptmp**2+E*log(temptmp))/100
indvap(it,iz,id)=indvap(it,iz,id)+1.


    endif
    enddo


A=Aw; B=Bw; C=Cw; D=Dw; E=Ew
    do it=32,tempmax  !! liq
    if( (temp(ix,iy,iz,id).ge.tempmod(it)).and.  &
        (temp(ix,iy,iz,id).lt.tempmod(it+1)).and.(temp(ix,iy,iz,id).ne.undef) &
       .and.(rhum(ix,iy,iz,id).gt.0.1) )then

temptmp=temp(ix,iy,iz,id)
vap(it,iz,id)=vap(it,iz,id)+rhum(ix,iy,iz,id)*  &
exp(A*(1/temptmp)+B+C*temptmp+D*temptmp**2+E*log(temptmp))/100
indvap(it,iz,id)=indvap(it,iz,id)+1.


    endif
    enddo



enddo
enddo
enddo
enddo

where(indvap>0.)
vap=vap/indvap
elsewhere
vap=-9999.
endwhere

print *, 'creation fichier'

call create_3d(trim(modfileout),'temp','pres','vap', &
               tempmax,altmax2,daymax,presnivs,tempmod,vap)
!subroutine create_3d(fname,ndim1,ndim2,ndim3,dim1,dim2,dim3,var1,var2,var3)


print *, 'end'

contains 

subroutine check(status)

implicit none
include 'netcdf.inc'   

 integer, intent ( in) :: status
    
    if(status /= nf90_noerr) then 
      print *, trim(nf90_strerror(status))
      stop "Stopped"
    end if
end subroutine check  

subroutine rdnc4(fname,lon,lat,alt,day,var,nvar)
  use netcdf
  implicit none
character(len=*)  ::  nvar, fname
integer  ::  lon,lat,alt,day
integer  ::  ncid,varid
real,dimension(lon,lat,alt,day)  :: var

call check(NF90_OPEN(trim(fname),NF90_NOWRITE,ncid))
call check(NF90_inq_varid(ncid,nvar,varid))
call check(NF90_get_var(ncid,varid,var))
call check(NF90_CLOSE(ncid))

end subroutine rdnc4


subroutine create_3d(fname,ndim1,ndim2,ndim3,dim1,dim2,dim3,var1,var2,var3)

    use netcdf
    implicit none

    character(LEN=*)   ::  fname,ndim1,ndim2,ndim3
    integer,parameter  ::  ndims=3
     real,parameter  ::  nan=-9999
 
    integer  ::  dim(ndims),nc_id,i,dim1,dim2,dim3
    integer  ::  varid, varid2,time_varid,varid3	 
    integer  ::  dimid1,dimid2,dimid3
   real,dimension(dim1)  ::  var1
   real,dimension(dim2)  ::  var2
   real*4,dimension(dim1,dim2,dim3)  ::  var3

   

! CREATION DU FICHIER nf90_create
    call check(nf90_create(fname, NF90_CLOBBER,  &
               nc_id))

! DEFINITION DES DIMENSIONS
    call check(nf90_def_dim(nc_id,ndim1, dim1, dimid1))
    call check(nf90_def_dim(nc_id,ndim2, dim2, dimid2 ))   
    call check(nf90_def_dim(nc_id,ndim3, dim3, dimid3))

! DIMENSION DES MATRICES  
    dim = (/dimid1, dimid2, dimid3/)

! VARIABLE 4D +time
    call check(nf90_def_var(nc_id, 'pres', NF90_FLOAT, dim2, varid))
    call check(nf90_put_att(nc_id, varid2, '_FillValue',nan))
    call check(nf90_def_var(nc_id, 'temp', NF90_FLOAT, dim1, varid2))
    call check(nf90_put_att(nc_id, varid, '_FillValue',nan))
    call check(nf90_def_var(nc_id, 'vap', NF90_FLOAT, dim, varid3))
    call check(nf90_put_att(nc_id, varid3, '_FillValue',nan))

! FIN DU MODE DEFINITION
    call check(nf90_enddef(nc_id))

! ENREGISTREMENT DES VARIABLES/VECTEURS DIMENSIONS
    call check(nf90_put_var(nc_id, varid, var1))
    call check(nf90_put_var(nc_id, varid2, var2))
    call check(nf90_put_var(nc_id, varid3, var3))

    call check(nf90_close(nc_id))

end subroutine create_3d

endprogram