       program ciesielski2couvreux

       IMPLICIT NONE

       INTEGER nbt, nlev
       INTEGER nbtf, nlevf
       INTEGER nbtb, nbtq, nlevq
       PARAMETER (nbt=736, nlev=40)
       PARAMETER (nbtf=736, nlevf=40)
       PARAMETER (nbtb=736)
       PARAMETER (nbtq=736, nlevq=40)
       INTEGER yy(nbt),mm(nbt),dd(nbt)
       INTEGER yyf(nbtf),mmf(nbtf),ddf(nbtf)
       INTEGER yyb(nbtb),mmb(nbtb),ddb(nbtb)
       INTEGER yyq(nbtq),mmq(nbtq),ddq(nbtq)
       INTEGER hh(nbt),hhf(nbtf),hhb(nbtb),hhq(nbtq)

       INTEGER k,l
       REAL p(nlev,nbt),z(nlev,nbt),u(nlev,nbt),v(nlev,nbt),omega(nlev,nbt)
       REAL T(nlev,nbt),theta(nlev,nbt),wmr(nlev,nbt),rh(nlev,nbt),div(nlev,nbt),vort(nlev,nbt)
       REAL hT(nlevf,nbtf), vT(nlevf,nbtf), hq(nlevf,nbtf), vq(nlevf,nbtf)
       REAL eo(nbtb),po2(nbtb),sh(nbtb),po1(nbtb),qr(nbtb),qrnet_ceres(nbtb)
       REAL ts(nbtb)

       REAL pf(nlevf,nbtf)
       REAL ug(nlevf,nbtf), vg(nlevf,nbtf)
       REAL dT(nlevf,nbtf), dq(nlevf,nbtf), dtrad(nlevf,nbtf)
       REAL dth(nlevf,nbtf), vth(nlevf,nbtf),hth(nlevf,nbtf)
       REAL dr(nlevf,nbtf), vr(nlevf,nbtf),hr(nlevf,nbtf)
       REAL du(nlevf,nbtf), vu(nlevf,nbtf),hu(nlevf,nbtf)
       REAL dv(nlevf,nbtf), vv(nlevf,nbtf),hv(nlevf,nbtf)

       REAL pq(nlevq,nbtq),q1(nlevq,nbtq),q2(nlevq,nbtq)
       REAL uw(nlevq,nbtq),vw(nlevq,nbtq)

!Lecture des champs moyens fournis par Ciesielski
       open(20,file='fields.ssa_3b')

       do l=1,nbt
          read (20,700) yy(l),mm(l),dd(l),hh(l)
          700 format (4(1x,i2))
          do k=1,nlev             
             read(20,800) p(k,l),z(k,l),u(k,l),v(k,l),omega(k,l),T(k,l),theta(k,l),wmr(k,l),rh(k,l),div(k,l),vort(k,l)
          800 format (11(1x,f8.2))
          enddo
       enddo

!Lecture des forçages grande-échelle fournis par Ciesielski
       open(21,file='lsf.ssa_3b')

       do l=1,nbtf
          read (21,700) yyf(l),mmf(l),ddf(l),hhf(l)
          do k=1,nlevf             
             read(21,900) pf(k,l),hT(k,l),vT(k,l),hq(k,l),vq(k,l)
             900 format (f8.2,1p,4e11.3) 
          enddo
       enddo

!Lecture des flux fournis par Ciesielski
       open(22,file='eopo.ssa_3b')

       do l=1,nbtb
          read(22,1000) yyb(l),mmb(l),ddb(l),hhb(l),eo(l),po2(l),sh(l),po1(l),qr(l),qrnet_ceres(l)
          1000 format (4i3,6f8.2)
       enddo

!Lecture des q1,q2 fournis par Ciesielski
       open(23,file='q1q2.nsa_3b')

        do l=1,nbtq
          read (23,700) yyq(l),mmq(l),ddq(l),hhq(l)
          do k=1,nlevq             
             read(23,1100) pq(k,l),q1(k,l),q2(k,l)
             1100 format (3f8.2) 
          enddo
       enddo

!on s'affranchit des undefined values (temporaire)
       do l=2,nbtf
           do k=1,nlevf
              if (hT(k,l).eq.-999) then
                 hT(k,l)=hT(k,l-1)
                 print*,'missing value hT',k,l
              endif
              if (vT(k,l).eq.-999) then
                 vT(k,l)=vT(k,l-1)
                 print*,'missing value vT',k,l
              endif
              if (omega(k,l).eq.-999) then
                 omega(k,l)=omega(k,l-1)
                 print*,'missing value w',k,l
              endif
           enddo
       enddo

!Initialisations des champs non fournis:
        do l=1,nbtf
           do k=1,nlevf
              ug(k,l)=0.
              vg(k,l)=0.
              hu(k,l)=0.
              vu(k,l)=0.
              hv(k,l)=0.
              vv(k,l)=0.
              uw(k,l)=0.
              vw(k,l)=0.
              dT(k,l)=hT(k,l)+vT(k,l)
              dtrad(k,l)=0.
              dq(k,l)=hq(k,l)+vq(k,l)
              du(k,l)=hu(k,l)+vu(k,l)
              dv(k,l)=hv(k,l)+vv(k,l)
            enddo
         enddo
         do l=1,nbtb
            ts(l)=T(1,l)+273.15
         enddo


!Ré-Ecriture au format LMDZ
       open(30,file='profil_moy.txt')
       write(30,*) nbt
       write(30,*) nlev
       do l=1,nbt
          write(30,*) yy(l),mm(l),dd(l),hh(l)*3600.
       do k=1,nlev
          write(30,*) p(k,l)*100.,z(k,l),u(k,l),v(k,l),T(k,l)+273.15,theta(k,l) &
     &               ,wmr(k,l)/1000./(1+wmr(k,l)/1000.),wmr(k,l)/1000.,rh(k,l)
       enddo
       enddo
       close(30)
     
       open(31,file='profil_adv.txt')
       write(31,*) "TEND",0,0,1,1,1
       write(31,*) "NUDG",7200,7200,0,0,0
       write(31,*) nbtf
       write(31,*) nlevf
       do l=1,nbtf
          write(31,*) yyf(l),mmf(l),ddf(l),hhf(l)*3600.
          do k=1,nlevf 
             write(31,*) pf(k,l)*100.,z(k,l),ug(k,l),vg(k,l),omega(k,l)*100./3600.,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)/1000.,hq(k,l)/1000.,vq(k,l)/1000.,dth(k,l),hth(k,l),vth(k,l) &
     &                   ,dr(k,l)/1000.,hr(k,l)/1000.,vr(k,l)/1000.
          enddo
       enddo
       close(31)

       open(33,file='uw_vw.txt')
       write(33,*) nbtf
       write(33,*) nlevf
       do l=1,nbtf
          write(33,*) yyf(l),mmf(l),ddf(l),hhf(l)*3600.
          do k=1,nlevf 
             write(33,*) pf(k,l)*100.,uw(k,l),vw(k,l)
          enddo
       enddo
       close(33)

       open(32,file='time_flux.txt') 
       write(32,*) nbtb
       do l=1,nbtb 
          write(32,*) yyb(l),mmb(l),ddb(l),hhb(l)*3600.,sh(l),eo(l),ts(l)
       enddo
       close(32)

       open(33,file='q1_q2.txt') 
       write(33,*) nbtq
       write(33,*) nlevq
       do l=1,nbtq
          write(33,*) yyq(l),mmq(l),ddq(l),hhq(l)*3600.
          do k=1,nlevq 
             write(33,*) pq(k,l)*100.,q1(k,l),q2(k,l)
          enddo
       enddo

       END 
