      program create_netcdf

      IMPLICIT none


      INTEGER llm,nbt
      INTEGER llh,nbth
      INTEGER llmf,nbtf
      INTEGER nbts
      INTEGER llmq,nbtq
      INTEGER llmw,nbtw
      INTEGER lm,nbstep
! llm, nbt: nb de niveaux et de pdt des champs moyens
! llmf, nbtf: nb de niveaux et de pdt des tendances
! nbts: nb de pdt des flux
! llmq, nbtq: nb de niveaux et de pdt de q1 et q2
! llmw, nbtw: nb de niveaux et de pdt de uw et vw
! lm,nbstep: dimensions du format netcdf

!     print*, 'Nombre de niveaux des champs moyens'
!     read(5,*) llm
!     print*, llm
!     print*, 'nombre de pas de temps des champs moyens'
!     read(5,*) nbt
!     print*, nbt
!     print*, 'Nombre de niveaux des advections'
!     read(5,*) llmf
!     print*, llmf
!     print*, 'nombre de pas de temps des advections'
!     read(5,*) nbtf
!     print*, nbtf
!     print*, 'nombre de pas de temps des flux'
!     read(5,*) nbts
!     print*, nbts
!     print*, 'Nombre de niveaux des q1 q2'
!     read(5,*) llmq
!     print*, llmq
!     print*, 'nombre de pas de temps des q1 q2'
!     read(5,*) nbtq
!     print*, nbtq
!
       llm=40
       nbt=1
       llh=41
       nbth=1
       llmf=40
       nbtf=31  ! forcing each 30 minutes during 15h30
       nbts=31
       llmq=40
       nbtq=31
       llmw=40
       nbtw=31
       lm=max(llm,llmf)
       print*,'avant prepare_armcu, lm=',lm
       nbstep=max(nbt,nbtf,nbts)
       print*,'avant prepare_armcu, nbstep=',nbstep
       CALL prepare_armcu(llm,nbt,llh,nbth,llmf,nbtf,nbts,llmq,nbtq,llmw,nbtw,lm,nbstep)

END


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      SUBROUTINE prepare_armcu(llm,nbt,llh,nbth,llmf,nbtf,nbts,llmq,nbtq,llmw,nbtw,lm,nbstep)
    

      IMPLICIT NONE

include "netcdf.inc"
include "YOMCST.h"
include "YOETHF.h"
      EXTERNAL suphel    ! initialiser certaines constantes

!Dimensions
!Fichiers etats moyens: llm, nbt
!Fichiers forçages: llmf, nbtf
!Fichiers flux de surface: nbts
!Fichiers de sortie: lm,nbstep

      INTEGER llm, nbt,nbt_autre
      INTEGER llmf, nbtf
      INTEGER nbts
      INTEGER llmq, nbtq
      INTEGER llmw, nbtw
      INTEGER lm,nbstep
      INTEGER l,k
      INTEGER llh,nbth

      INTEGER xm,ym
      PARAMETER(xm=1,ym=1)
!Variables lues
!Dates
      REAL*8 year(nbt),month(nbt),day(nbt)
      REAL*8 hour(nbt),time(nbt)
!     REAL*8 yearf(nbtf),monthf(nbtf),dayf(nbtf)
!     REAL*8 hourf(nbtf),timef(nbtf)
!     REAL*8 years(nbts),months(nbts),days(nbts)
!     REAL*8 hours(nbts),times(nbts)
!     REAL*8 yearq(nbtq),monthq(nbtq),dayq(nbtq)
!     REAL*8 hourq(nbtq),timeq(nbtq)
      REAL*8 time_out(nbstep),poub
      REAL*8 time_out_min,time_out_max
      
!Définition lon, lat
      REAL*8 lon(xm),lat(ym),zrho,missing_val

!niveaux verticaux
      REAL*8 zz(llm,nbt),zzf(llmf,nbtf),zlev(lm,nbstep)
      REAL*8 ap(llh),bp(llh),pph(llh),zzh(llh)

!profils moyens lus
      REAL*8 th(llm,nbt),mwr(llm,nbt),u(llm,nbt),v(llm,nbt)
      REAL*8 pp(llm,nbt),qv(llm,nbt),ql(llm,nbt),qi(llm,nbt),T(llm,nbt),rh(llm,nbt)
      REAL*8 vent,uw(llmw,nbtw),vw(llmw,nbtw),ppw(llmw,nbtw)
      

!forcages grande echelle lus
      REAL*8 ppf(llmf,nbtf),rho(llmf,nbtf)
      REAL*8 ug(llmf,nbtf),vg(llmf,nbtf),w(llmf,nbtf),omega(llmf,nbtf)
      REAL*8 du(llmf,nbtf),hu(llmf,nbtf),vu(llmf,nbtf)
      REAL*8 dv(llmf,nbtf),hv(llmf,nbtf),vv(llmf,nbtf)
      REAL*8 dT(llmf,nbtf),hT(llmf,nbtf),vT(llmf,nbtf),dtrad(llmf,nbtf)
      REAL*8 dq(llmf,nbtf),hq(llmf,nbtf),vq(llmf,nbtf)
      REAL*8 dth(llmf,nbtf),hth(llmf,nbtf),vth(llmf,nbtf)
      REAL*8 dr(llmf,nbtf),hr(llmf,nbtf),vr(llmf,nbtf)

!q1 q2 lus
      REAL*8 ppq(llmq,nbtq)
      REAL*8 qq1(llmq,nbtq),qq2(llmq,nbtq)

!Paramètres de forçages
      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)

      REAL*8 ustar

!Variables intermédaires
      REAL*8 ppf_tmpv(lm,nbtf)
      REAL*8 ug_tmpv(lm,nbtf),vg_tmpv(lm,nbtf),w_tmpv(lm,nbtf)
      REAL*8 du_tmpv(lm,nbtf),hu_tmpv(lm,nbtf),vu_tmpv(lm,nbtf)
      REAL*8 dv_tmpv(lm,nbtf),hv_tmpv(lm,nbtf),vv_tmpv(lm,nbtf)
      REAL*8 dT_tmpv(lm,nbtf),hT_tmpv(lm,nbtf),vT_tmpv(lm,nbtf)
      REAL*8 dtrad_tmpv(lm,nbtf)
      REAL*8 dq_tmpv(lm,nbtf),hq_tmpv(lm,nbtf),vq_tmpv(lm,nbtf)
      REAL*8 dth_tmpv(lm,nbtf),hth_tmpv(lm,nbtf),vth_tmpv(lm,nbtf)
      REAL*8 dr_tmpv(lm,nbtf),hr_tmpv(lm,nbtf),vr_tmpv(lm,nbtf)      

      REAL*8 th_tmpv(lm,nbt),mwr_tmpv(lm,nbt),u_tmpv(lm,nbt),v_tmpv(lm,nbt)
      REAL*8 pp_tmpv(lm,nbt),qv_tmpv(lm,nbt),ql_tmpv(lm,nbt),qi_tmpv(lm,nbt),T_tmpv(lm,nbt),rh_tmpv(lm,nbt)

      REAL*8 ppq_tmpv(lm,nbtq),q1_tmpv(lm,nbtq),q2_tmpv(lm,nbtq)

      REAL*8 zz_tmpv(lm,nbt),zzf_tmpv(lm,nbtf)

!Variables ecrites dans fichier netcdf aux dimensions lm,nbstep
!profils moyen
      REAL*8 vitu(lm,nbstep),vitv(lm,nbstep)
      REAL*8 pres(lm,nbstep),ovap(lm,nbstep),temp(lm,nbstep),rhum(lm,nbstep)
!     REAL*8 theta(lm,nbstep),rv(lm,nbstep)     
      REAL*8 theta(lm,nbstep)     

!flux de surface
      REAL*8 flat(nbstep),sens(nbstep),tsurf(nbstep)

!flux de surface lus
      REAL*8 flxs(lm,nbstep),flxl(lm,nbstep),ts(nbstep),ps(nbstep),orog

!forcages grande echelle
      REAL*8 vitug(lm,nbstep),vitvg(lm,nbstep),vitw(lm,nbstep)
      REAL*8 advT(lm,nbstep),T_advh(lm,nbstep),T_advv(lm,nbstep),radT(lm,nbstep)
      REAL*8 advq(lm,nbstep),q_advh(lm,nbstep),q_advv(lm,nbstep)   
      REAL*8 advth(lm,nbstep),th_advh(lm,nbstep),th_advv(lm,nbstep)
      REAL*8 advr(lm,nbstep),r_advh(lm,nbstep),r_advv(lm,nbstep)   
      REAL*8 advu(lm,nbstep),u_advh(lm,nbstep),u_advv(lm,nbstep)
      REAL*8 advv(lm,nbstep),v_advh(lm,nbstep),v_advv(lm,nbstep)  

      REAL*8 q1(lm,nbstep),q2(lm,nbstep)  

      REAL*8 timestep,dtm,dtf,dts
      REAL*8 day_init,month_init
      REAL*8 dayf_init,monthf_init
      REAL*8 days_init,months_init
      INTEGER nday,ndayf,ndays,ndayq

!Définition des sorties
      integer nbvar3d_out
      parameter (nbvar3d_out=39)
      character (len=50), dimension(nbvar3d_out) :: varname3d_out

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

      integer var3didout(39),toto(lm)

!     character*4 poub
      real*8 ppp(lm)

      REAL*8 pzero
      PARAMETER(pzero=1.e5)

      CALL suphel
      print *,'Apres CALL suphel'

!     Read in data: profil.txt already exists in
!     http://www.knmi.nl/samenw/eurocs/ARM/case_ARM_html/case_ARM.html

      missing_val=-9999
      ts(1)=296.4
      ps(1)=97000.
      orog=0.
      ustar=missing_val

      open(10,file='profil.txt')
      open(20,file='profil_moy.txt')
      read(10,*) nbt
      write(20,*) nbtf
      read(10,*) llm
      write(20,*) llm
      read (10,*) year(1),month(1),day(1),hour(1)
      write (20,*) year(1),month(1),day(1),hour(1)
!     print *,'year=',year(1),month(1),day(1),hour(1)
! Attention !! rh est en fait rv
      do k=1,llm
         read(10,*) zz(k,1),pp(k,1),u(k,1),v(k,1),th(k,1),T(k,1),qv(k,1),rh(k,1)
         pp(k,1)=pp(k,1)*100.
         qv(k,1)=qv(k,1)/1000.
         ql(k,1)=missing_val
         qi(k,1)=missing_val
         rh(k,1)=rh(k,1)/1000.
         write(20,*) pp(k,1),zz(k,1),u(k,1),v(k,1),T(k,1),th(k,1),qv(k,1),ql(k,1),qi(k,1),rh(k,1)
         print *,'data=',pp(k,1),zz(k,1),T(k,1),th(k,1),qv(k,1)
      enddo
! Meme s il n y a qu un seul profil initial, on le repete autant que les
! forcages (nbtf)
      do l=2,nbtf
          write(20,*) year(1),month(1),day(1),hour(1)+(l-1)*1800.
          do k=1,llm
           write(20,*) pp(k,1),zz(k,1),u(k,1),v(k,1),T(k,1),th(k,1),qv(k,1),ql(k,1),qi(k,1),rh(k,1)
          enddo
      enddo
      close(10)
      close(20)
      print *,'Profil initial: nbt,llm',nbt,llm

!     Read in data: cab.txt already exists in
!     http://www.knmi.nl/samenw/eurocs/ARM/case_ARM_html/case_ARM.html

      open(10,file='cab.txt')
      open(20,file='coord_a_b.txt')
      read(10,*) nbth
      write(20,*) nbth
      read(10,*) llh
      write(20,*) llh
      do k=1,llh
         read(10,*) l,ap(k),bp(k)
         pph(k)= ap(k)+bp(k)*ps(1)
         write(20,*) ap(k),bp(k),pph(k)
!        print *,'ap,bp,pph=',k,ap(k),bp(k),pph(k)
      enddo
      close(10)
      close(20)
      print *,'Profil initial: nbth,llh',nbth,llh
!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(10,file='forc_orig.txt')
      open(20,file='profil_adv.txt')
! Combien de profils de forcing : toutes les demi heures pendant 24h
      read(10,*) nbtf
      write(20,*) nbtf
! Combien de niveaux : autant que pour les profils initiaux
      write(20,*) llmf
!
      do l=1,nbtf
      do k=1,llmf
         ppf(k,l)=pp(k,1)
         zzf(k,l)=zz(k,1)
!        print *,'pp,zz=',k,l,ppf(k,l),zzf(k,l),zz(k,1)
      enddo
      enddo
! Lecture des profils
      do l=1,nbtf
          read(10,*) year(l),month(l),day(l),hour(l),flxs(1,l),flxl(1,l),advT(1,l),radT(1,l),advq(1,l)
          write(20,*) year(l),month(l),day(l),hour(l)
          do k=1,llmf
!The zonal u-component of the geostrophic wind is decreasing with 2.0 * 10^-3s-1 corresponding with the observed wind
             ug(k,l)=u(k,1)
!The geostrophic v-component is assumed to be equal to the meridional wind v. 
             vg(k,l)=v(k,1)
! w en m/s
             w(k,l)=0.
             rho(k,l)=ppf(k,l)/T(k,1)/RD    ! conversion m/s en Pa/s
             omega(k,l)=-1.*w(k,l)*rho(k,l)*RG
             du(k,l)=missing_val
             hu(k,l)=missing_val
             vu(k,l)=missing_val
             dv(k,l)=missing_val
             hv(k,l)=missing_val
             vv(k,l)=missing_val
!            dT(k,l)=missing_val
             hT(k,l)=missing_val
             vT(k,l)=missing_val
! Large Scale Temperature Forcing due to advection and radiation, dT/dt [K/h]
! and dq/dt [g/kg/h]
             if (zzf(k,l).LT.1000.) THEN 
                dth(k,l)= (advT(1,l)+radT(1,l))/3600.
                dq(k,l)= advq(1,l)/1000./3600.
             elseif ((zzf(k,l).GE.1000.).AND.(zzf(k,l).LE.3000.)) THEN 
                dth(k,l)= (advT(1,l)+radT(1,l))/3600.
                dth(k,l)=dth(k,l)*(1.-(zzf(k,l)-1000.)/2000.)
                dq(k,l)= advq(1,l)/1000./3600.
                dq(k,l)=dq(k,l)*(1.-(zzf(k,l)-1000.)/2000.)
             else 
                dth(k,l)=0.
                dq(k,l)=0.
             endif
!----------------------------------------------
! Le cas est force est theta, on convertit dth,hth,vth
! (adv totale, adv horizontale et adv verticale) en dt,ht,vt
! teta=temp*(pzero/play)**rkappa
!----------------------------------------------
             dT(k,l)=dth(k,l)/(pzero/ppf(k,l))**rkappa
!Large Scale Horizontal Moisture Advection dqv/dt [(g/kg)/s] 
!               dq(k,l)=0.
                hq(k,l)=missing_val
                vq(k,l)=missing_val
                vt(k,l)=missing_val
!               dth(k,l)=0.
                hth(k,l)=missing_val
                vth(k,l)=missing_val
                dr(k,l)=missing_val
                hr(k,l)=missing_val
                vr(k,l)=missing_val
                dtrad(k,l)=missing_val
             write(20,*) ppf(k,l),zzf(k,l),ug(k,l),vg(k,l),w(k,l),omega(k,l),du(k,l),hu(k,l),vu(k,l) &
     &                 ,dv(k,l),hv(k,l),vv(k,l),dT(k,l),hT(k,l),vT(k,l),dtrad(k,l),dq(k,l) &
     &                 ,hq(k,l),vq(k,l),dth(k,l),hth(k,l),vth(k,l),dr(k,l),hr(k,l),vr(k,l) 
          enddo
      enddo
      close(10)
      close(20)

!Lecture des flux
      open(10,file='forc_orig.txt')
      open(20,file='time_flux.txt')
! Combien de flux?
      read(10,*) nbtf
      write(20,*) nbtf
      do l=1,nbtf
         read(10,*) year(l),month(l),day(l),hour(l),flxs(1,l),flxl(1,l)
         zrho= pp(1,1)/rd/T(1,1)
         write(20,*) year(l),month(l),day(l),hour(l),flxs(1,l),flxl(1,l),ts(1),ustar,ps(1),orog
!        print *,'pp,T,sens,flat',pp(1,l),T(1,l),flxs(1,l),flxl(1,l)
      enddo
      close(10)
      close(20)

!     read Q1, Q2

      open(20,file='q1_q2.txt')
! Combien de profils de q1 q2 ?
      write(20,*),nbtq
! Combien de niveaux ?
      write(20,*),llmq
      do l=1,nbtq
          write (20,*) year(l),month(l),day(l),hour(l)
          do k=1,llmq
             ppq(k,l)=pp(k,1)
             qq1(k,l)=missing_val
             qq2(k,l)=missing_val
             write(20,*) ppq(k,l),qq1(k,l),qq2(k,l)
          enddo
      enddo
      close(20)

      open(20,file='uw_vw.txt')
! Combien de profils de q1 q2 ?
      write(20,*),nbtq
! Combien de niveaux ?
      write(20,*),llmq
      do l=1,nbtq
          write (20,*) year(l),month(l),day(l),hour(l)
          do k=1,llmq
             ppw(k,l)=pp(k,1)
             uw(k,l)=missing_val
             vw(k,l)=missing_val
             write(20,*) ppw(k,l),uw(k,l),vw(k,l)
          enddo
      enddo
      close(20)

     stop
     end
