!----------------------------------------------------------------------
! forcing_les = .T. : Impose a constant cooling 
! forcing_radconv = .T. : Pure radiative-convective equilibrium:
!----------------------------------------------------------------------

      if (forcing_les .or. forcing_radconv .or. forcing_GCSSold ) then

!----------------------------------------------------------------------
! Read profiles from files: prof.inp.001 and lscale.inp.001
! (repris de readlesfiles)
!----------------------------------------------------------------------

      call readprofiles(nlev_max,kmax,height,
     .           tttprof,qtprof,uprof,vprof,
     .           e12prof,ugprof,vgprof,
     .           wfls,dqtdxls,dqtdyls,dqtdtls,
     .           thlpcar) 

! compute altitudes of play levels.
      zlay(1) =zsurf +  rd*tsurf*(psurf-play(1))/(rg*psurf)
      do l = 2,llm
        zlay(l) = zlay(l-1)+rd*tsurf*(psurf-play(1))/(rg*psurf)
      enddo

!----------------------------------------------------------------------
! Interpolation of the profiles given on the input file to
! model levels
!----------------------------------------------------------------------
      zlay(1) = zsurf +  rd*tsurf*(psurf-play(1))/(rg*psurf)
      do l=1,llm
        ! Above the max altutide of the input file

        if (zlay(l)<height(kmax)) mxcalc=l

        frac = (height(kmax)-zlay(l))/(height (kmax)-height(kmax-1))
        ttt =tttprof(kmax)-frac*(tttprof(kmax)-tttprof(kmax-1))
        if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile
          temp(l) = ttt*(play(l)/pzero)**rkappa
          teta(l) = ttt
        else
          temp(l) = ttt
          teta(l) = ttt*(pzero/play(l))**rkappa
        endif
          print *,' temp,teta ',l,temp(l),teta(l)
        q(l,1)  = qtprof(kmax)-frac*( qtprof(kmax)- qtprof(kmax-1))
        u(l)    =  uprof(kmax)-frac*(  uprof(kmax)-  uprof(kmax-1))
        v(l)    =  vprof(kmax)-frac*(  vprof(kmax)-  vprof(kmax-1))
        ug(l)   = ugprof(kmax)-frac*( ugprof(kmax)- ugprof(kmax-1))
        vg(l)   = vgprof(kmax)-frac*( vgprof(kmax)- vgprof(kmax-1))
        omega(l)=   wfls(kmax)-frac*(   wfls(kmax)-   wfls(kmax-1))

        dq_dyn(l,1) = dqtdtls(kmax)-frac*(dqtdtls(kmax)-dqtdtls(kmax-1))
        dt_cooling(l)
     .          =thlpcar(kmax)-frac*(thlpcar(kmax)-thlpcar(kmax-1))
        do k=2,kmax
          frac = (height(k)-zlay(l))/(height(k)-height(k-1))
          if(l==1) print*,'k, height, tttprof',k,height(k),tttprof(k)
          if(zlay(l)>height(k-1).and.zlay(l)<height(k)) then
            ttt =tttprof(k)-frac*(tttprof(k)-tttprof(k-1))
        if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile
          temp(l) = ttt*(play(l)/pzero)**rkappa
          teta(l) = ttt
        else
          temp(l) = ttt
          teta(l) = ttt*(pzero/play(l))**rkappa
        endif
          print *,' temp,teta ',l,temp(l),teta(l)
            q(l,1)  = qtprof(k)-frac*( qtprof(k)- qtprof(k-1))
            u(l)    =  uprof(k)-frac*(  uprof(k)-  uprof(k-1))
            v(l)    =  vprof(k)-frac*(  vprof(k)-  vprof(k-1))
            ug(l)   = ugprof(k)-frac*( ugprof(k)- ugprof(k-1))
            vg(l)   = vgprof(k)-frac*( vgprof(k)- vgprof(k-1))
            omega(l)=   wfls(k)-frac*(   wfls(k)-   wfls(k-1))
            dq_dyn(l,1)=dqtdtls(k)-frac*(dqtdtls(k)-dqtdtls(k-1))
            dt_cooling(l)
     .              =thlpcar(k)-frac*(thlpcar(k)-thlpcar(k-1))
          elseif(zlay(l)<height(1)) then ! profils uniformes pour z<height(1)
            ttt =tttprof(1)
        if (forcing_GCSSold .AND. tp_ini_GCSSold) then ! pot. temp. in initial profile
          temp(l) = ttt*(play(l)/pzero)**rkappa
          teta(l) = ttt
        else
          temp(l) = ttt
          teta(l) = ttt*(pzero/play(l))**rkappa
        endif
            q(l,1)  = qtprof(1)
            u(l)    =  uprof(1)
            v(l)    =  vprof(1)
            ug(l)   = ugprof(1)
            vg(l)   = vgprof(1)
            omega(l)=   wfls(1)
            dq_dyn(l,1)  =dqtdtls(1)
            dt_cooling(l)=thlpcar(1)
          endif
        enddo 

        temp(l)=max(min(temp(l),350.),150.)
        rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
        if (l .lt. llm) then
          zlay(l+1) = zlay(l) + (play(l)-play(l+1))/(rg*rho(l))
        endif
        omega2(l)=-rho(l)*omega(l)
        omega(l)= omega(l)*(-rg*rho(l)) !en Pa/s
        if (l>1) then
        if(zlay(l-1)>height(kmax)) then
           omega(l)=0.0
           omega2(l)=0.0
        endif   
        endif
        if(q(l,1)<0.) q(l,1)=0.0
        q(l,2)  = 0.0
      enddo

      endif ! forcing_les .or. forcing_GCSSold
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------
! Forcing for GCSSold:
!---------------------------------------------------------------------
      if (forcing_GCSSold) then
       fich_gcssold_ctl = './forcing.ctl'
       fich_gcssold_dat = './forcing8.dat'
       call copie(llm,play,psurf,fich_gcssold_ctl)
       call get_uvd2(it,timestep,fich_gcssold_ctl,fich_gcssold_dat,
     :               ht_gcssold,hq_gcssold,hw_gcssold,
     :               hu_gcssold,hv_gcssold,
     :               hthturb_gcssold,hqturb_gcssold,Ts_gcssold,
     :               imp_fcg_gcssold,ts_fcg_gcssold,
     :               Tp_fcg_gcssold,Turb_fcg_gcssold)
       print *,' get_uvd2 -> hqturb_gcssold ',hqturb_gcssold
      endif ! forcing_GCSSold
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------
! Forcing for RICO:
!---------------------------------------------------------------------
      if (forcing_rico) then

!       call writefield_phy('omega', omega,llm+1)
      fich_rico = 'rico.txt'
       call read_rico(fich_rico,nlev_rico,ps_rico,play
     :             ,ts_rico,t_rico,q_rico,u_rico,v_rico,w_rico
     :             ,dth_rico,dqh_rico)
        print*, ' on a lu et prepare RICO'

       mxcalc=llm
       print *, airefi, ' airefi '
       do l = 1, llm
       rho(l)  = play(l)/(rd*t_rico(l)*(1.+(rv/rd-1.)*q_rico(l)))
       temp(l) = t_rico(l)
       q(l,1) = q_rico(l)
       q(l,2) = 0.0
       u(l) = u_rico(l)
       v(l) = v_rico(l)
       ug(l)=u_rico(l)
       vg(l)=v_rico(l)
       omega(l) = -w_rico(l)*rg
       omega2(l) = omega(l)/rg*airefi
       enddo
      endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------
! Forcing from TOGA-COARE experiment (Ciesielski et al. 2002) :
!---------------------------------------------------------------------

      if (forcing_toga) then

! read TOGA-COARE forcing (native vertical grid, nt_toga timesteps):
      fich_toga = './d_toga/ifa_toga_coare_v21_dime.txt'
      CALL read_togacoare(fich_toga,nlev_toga,nt_toga
     :         ,ts_toga,plev_toga,t_toga,q_toga,u_toga,v_toga,w_toga
     :         ,ht_toga,vt_toga,hq_toga,vq_toga)

       write(*,*) 'Forcing TOGA lu'

! time interpolation for initial conditions:
      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
      CALL interp_toga_time(daytime,day1,annee_ref
     i             ,year_ini_toga,day_ju_ini_toga,nt_toga,dt_toga
     i             ,nlev_toga,ts_toga,plev_toga,t_toga,q_toga,u_toga
     i             ,v_toga,w_toga,ht_toga,vt_toga,hq_toga,vq_toga
     o             ,ts_prof,plev_prof,t_prof,q_prof,u_prof,v_prof,w_prof
     o             ,ht_prof,vt_prof,hq_prof,vq_prof)

! vertical interpolation:
      CALL interp_toga_vertical(play,nlev_toga,plev_prof
     :         ,t_prof,q_prof,u_prof,v_prof,w_prof
     :         ,ht_prof,vt_prof,hq_prof,vq_prof
     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
       write(*,*) 'Profil initial forcing TOGA interpole'

! initial and boundary conditions :
      tsurf = ts_prof
      write(*,*) 'SST initiale: ',tsurf
      do l = 1, llm
       temp(l) = t_mod(l)
       q(l,1) = q_mod(l)
       q(l,2) = 0.0
       u(l) = u_mod(l)
       v(l) = v_mod(l)
       omega(l) = w_mod(l)
       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq
!?       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
!?       omega2(l)=-rho(l)*omega(l)
       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
       d_th_adv(l) = alpha*omega(l)/rcpd-(ht_mod(l)+vt_mod(l))
       d_q_adv(l,1) = -(hq_mod(l)+vq_mod(l))
       d_q_adv(l,2) = 0.0
      enddo

      endif ! forcing_toga
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------
! Forcing from TWPICE experiment (Shaocheng et al. 2010) :
!---------------------------------------------------------------------

      if (forcing_twpice) then
!read TWP-ICE forcings
      fich_twpice=
     : 'd_twpi/twp180iopsndgvarana_v2.1_C3.c1.20060117.000000.cdf'
      call read_twpice(fich_twpice,nlev_twpi,nt_twpi
     :     ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi
     :     ,ht_twpi,vt_twpi,hq_twpi,vq_twpi)

      write(*,*) 'Forcing TWP-ICE lu'
!Time interpolation for initial conditions using TOGA interpolation routine
         write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',daytime,day1   
      CALL interp_toga_time(daytime,day1,annee_ref
     i          ,year_ini_twpi,day_ju_ini_twpi,nt_twpi,dt_twpi,nlev_twpi
     i             ,ts_twpi,plev_twpi,t_twpi,q_twpi,u_twpi,v_twpi,w_twpi
     i             ,ht_twpi,vt_twpi,hq_twpi,vq_twpi
     o             ,ts_proftwp,plev_proftwp,t_proftwp,q_proftwp
     o             ,u_proftwp,v_proftwp,w_proftwp
     o             ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp)

! vertical interpolation using TOGA interpolation routine:
!      write(*,*)'avant interp vert', t_proftwp
      CALL interp_toga_vertical(play,nlev_twpi,plev_proftwp
     :         ,t_proftwp,q_proftwp,u_proftwp,v_proftwp,w_proftwp
     :         ,ht_proftwp,vt_proftwp,hq_proftwp,vq_proftwp
     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
!       write(*,*) 'Profil initial forcing TWP-ICE interpole',t_mod

! initial and boundary conditions :
!      tsurf = ts_proftwp
      write(*,*) 'SST initiale: ',tsurf
      do l = 1, llm
       temp(l) = t_mod(l)
       q(l,1) = q_mod(l)
       q(l,2) = 0.0
       u(l) = u_mod(l)
       v(l) = v_mod(l)
       omega(l) = w_mod(l)
       omega2(l)=omega(l)/rg*airefi ! flxmass_w calcule comme ds physiq

       alpha = rd*temp(l)*(1.+(rv/rd-1.)*q(l,1))/play(l)
!on applique le forcage total au premier pas de temps
!attention: signe different de toga
       d_th_adv(l) = alpha*omega(l)/rcpd+(ht_mod(l)+vt_mod(l))
       d_q_adv(l,1) = (hq_mod(l)+vq_mod(l))
       d_q_adv(l,2) = 0.0
      enddo     
       
      endif !forcing_twpice
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!---------------------------------------------------------------------
! Forcing from Arm_Cu case                   
! For this case, ifa_armcu.txt contains sensible, latent heat fluxes
! large scale advective forcing,radiative forcing 
! and advective tendency of theta and qt to be applied 
!---------------------------------------------------------------------

      if (forcing_armcu) then
! read armcu forcing :
       write(*,*) 'Avant lecture Forcing Arm_Cu'
      fich_armcu = './ifa_armcu.txt'
      CALL read_armcu(fich_armcu,nlev_armcu,nt_armcu,
     : sens_armcu,flat_armcu,adv_theta_armcu,
     : rad_theta_armcu,adv_qt_armcu)
       write(*,*) 'Forcing Arm_Cu lu'

!----------------------------------------------------------------------
! Read profiles from file: prof.inp.19 or prof.inp.40 
! For this case, profiles are given for two vertical resolution
! 19 or 40 levels
!
! Comment from: http://www.knmi.nl/samenw/eurocs/ARM/profiles.html
! Note that the initial profiles contain no liquid water! 
! (so potential temperature can be interpreted as liquid water 
! potential temperature and water vapor as total water)
! profiles are given at full levels
!----------------------------------------------------------------------

      call readprofile_armcu(nlev_max,kmax,height,play_mod,u_mod,
     .           v_mod,theta_mod,t_mod,qv_mod,rv_mod,ap,bp)

! time interpolation for initial conditions:
      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1

      print *,'Avant interp_armcu_time'
      print *,'daytime=',daytime
      print *,'day1=',day1
      print *,'annee_ref=',annee_ref
      print *,'year_ini_armcu=',year_ini_armcu
      print *,'day_ju_ini_armcu=',day_ju_ini_armcu
      print *,'nt_armcu=',nt_armcu
      print *,'dt_armcu=',dt_armcu
      print *,'nlev_armcu=',nlev_armcu
      CALL interp_armcu_time(daytime,day1,annee_ref
     i            ,year_ini_armcu,day_ju_ini_armcu,nt_armcu,dt_armcu
     i            ,nlev_armcu,sens_armcu,flat_armcu,adv_theta_armcu
     i            ,rad_theta_armcu,adv_qt_armcu,sens_prof,flat_prof
     i            ,adv_theta_prof,rad_theta_prof,adv_qt_prof)
       write(*,*) 'Forcages interpoles dans temps'

! No vertical interpolation if nlev imposed to 19 or 40
! The vertical grid stops at 4000m # 600hPa
      mxcalc=llm

! initial and boundary conditions :
!     tsurf = ts_prof
! tsurf read in lmdz1d.def
      write(*,*) 'Tsurf initiale: ',tsurf
      do l = 1, llm
       play(l)=play_mod(l)*100.
       presnivs(l)=play(l)
       zlay(l)=height(l)
       temp(l) = t_mod(l)
       teta(l)=theta_mod(l)
       q(l,1) = qv_mod(l)/1000.
! No liquid water in the initial profil
       q(l,2) = 0.
       u(l) = u_mod(l)
       ug(l)= u_mod(l)
       v(l) = v_mod(l)
       vg(l)= v_mod(l)
! Advective forcings are given in K or g/kg ... per HOUR
!      IF(height(l).LT.1000) THEN
!        d_th_adv(l) = (adv_theta_prof + rad_theta_prof)/3600.
!        d_q_adv(l,1) = adv_qt_prof/1000./3600.
!        d_q_adv(l,2) = 0.0
!      ELSEIF (height(l).GE.1000.AND.height(l).LT.3000) THEN
!        d_th_adv(l) = (adv_theta_prof + rad_theta_prof)*
!    :               (1-(height(l)-1000.)/2000.)
!        d_th_adv(l) = d_th_adv(l)/3600.
!        d_q_adv(l,1) = adv_qt_prof*(1-(height(l)-1000.)/2000.)
!        d_q_adv(l,1) = d_q_adv(l,1)/1000./3600.
!        d_q_adv(l,2) = 0.0
!      ELSE
!        d_th_adv(l) = 0.0
!        d_q_adv(l,1) = 0.0
!        d_q_adv(l,2) = 0.0
!      ENDIF
      enddo
! plev at half levels is given in proh.inp.19 or proh.inp.40 files
      plev(1)= ap(llm+1)+bp(llm+1)*psurf
      do l = 1, llm
      plev(l+1) = ap(llm-l+1)+bp(llm-l+1)*psurf
      print *,'Read_forc: l height play plev zlay temp',
     :   l,height(l),play(l),plev(l),zlay(l),temp(l)
      enddo
! For this case, fluxes are imposed
       fsens=-1*sens_prof
       flat=-1*flat_prof

      endif ! forcing_armcu
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!