Index: LMDZ5/branches/testing/libf/phylmd/1DUTILS.h
===================================================================
--- LMDZ5/branches/testing/libf/phylmd/1DUTILS.h	(revision 2056)
+++ LMDZ5/branches/testing/libf/phylmd/1DUTILS.h	(revision 2160)
@@ -1,3 +1,3 @@
-#include "../dyn3d/conf_gcm.F"
+#include "../dyn3d/conf_gcm.F90"
 #include "../dyn3d_common/q_sat.F"
 
@@ -140,4 +140,11 @@
        CALL getin('ok_flux_surf',ok_flux_surf)
 
+!Config  Key  = ok_prescr_ust
+!Config  Desc = ustar impose ou non
+!Config  Def  = false
+!Config  Help = ustar impose ou non
+       ok_prescr_ust = .false.
+       CALL getin('ok_prescr_ust',ok_prescr_ust)
+
 !Config  Key  = ok_old_disvert
 !Config  Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
@@ -355,5 +362,5 @@
       CHARACTER*(*) fichnom
 !Al1 plev tronque pour .nc mais plev(klev+1):=0
-      real :: plev(klon,klev),play (klon,klev),phi(klon,klev)
+      real :: plev(klon,klev+1),play (klon,klev),phi(klon,klev)
       real :: presnivs(klon,klev)
       real :: ucov(klon,klev),vcov(klon,klev),temp(klon,klev)
@@ -438,7 +445,4 @@
 !  Lecture des champs
 !
-      plev(1,klev+1)=0.
-      CALL get_field("plev",plev,found)
-      IF (.NOT. found) PRINT*, modname//'Le champ <Plev> est absent'
       CALL get_field("play",play,found)
       IF (.NOT. found) PRINT*, modname//'Le champ <Play> est absent'
@@ -457,4 +461,7 @@
       CALL get_field("omega2",omega2,found)
       IF (.NOT. found) PRINT*, modname//'Le champ <omega2> est absent'
+      plev(1,klev+1)=0.
+      CALL get_field("plev",plev(:,1:klev),found)
+      IF (.NOT. found) PRINT*, modname//'Le champ <Plev> est absent'
 
       Do iq=1,nqtot
@@ -822,5 +829,4 @@
       INTEGER np,ierr
       REAL pi,x
-      REAL :: scaleheight=8.
  
 !-----------------------------------------------------------------------
@@ -925,14 +931,8 @@
       PRINT *,  ap
  
-   print*,'scaleheight=',scaleheight
-  DO l = 1, llm
-     dpres(l) = bp(l) - bp(l+1)
-     presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
-     write(*, *)'PRESNIVS(', l, ')=', presnivs(l), ' Z ~ ', &
-          log(preff/presnivs(l))*scaleheight &
-          , ' DZ ~ ', scaleheight*log((ap(l)+bp(l)*preff)/ &
-          max(ap(l+1)+bp(l+1)*preff, 1.e-10))
-  ENDDO
-
+      DO l = 1, llm
+       dpres(l) = bp(l) - bp(l+1)
+       presnivs(l) = 0.5 *( ap(l)+bp(l)*preff + ap(l+1)+bp(l+1)*preff )
+      ENDDO
  
       PRINT *,' PRESNIVS '
@@ -2108,5 +2108,5 @@
         elseif(4000 < zlay(k) .and. zlay(k) < 9000) then
           q_rico(k)=1.8 + (0 - 1.8) /                                      &
-     &             (10000 - 4000) * (zlay(k) - 4000)
+     &             (9000 - 4000) * (zlay(k) - 4000)
         else
           q_rico(k)=0.0
@@ -2420,4 +2420,127 @@
           end
  
+!=====================================================================
+       SUBROUTINE interp_dice_vertical(play,nlev_dice,nt_dice,plev_prof   &
+     &         ,th_prof,qv_prof,u_prof,v_prof,o3_prof                     &
+     &         ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof         &
+     &         ,th_mod,qv_mod,u_mod,v_mod,o3_mod                          &
+     &         ,ht_mod,hq_mod,hu_mod,hv_mod,w_mod,omega_mod,mxcalc)
+ 
+       implicit none
+ 
+#include "dimensions.h"
+
+!-------------------------------------------------------------------------
+! Vertical interpolation of Dice forcing data onto model levels
+!-------------------------------------------------------------------------
+ 
+       integer nlevmax
+       parameter (nlevmax=41)
+       integer nlev_dice,mxcalc,nt_dice
+ 
+       real play(llm), plev_prof(nlev_dice) 
+       real th_prof(nlev_dice),qv_prof(nlev_dice)
+       real u_prof(nlev_dice),v_prof(nlev_dice) 
+       real o3_prof(nlev_dice)
+       real ht_prof(nlev_dice),hq_prof(nlev_dice)
+       real hu_prof(nlev_dice),hv_prof(nlev_dice)
+       real w_prof(nlev_dice),omega_prof(nlev_dice)
+ 
+       real th_mod(llm),qv_mod(llm)
+       real u_mod(llm),v_mod(llm), o3_mod(llm)
+       real ht_mod(llm),hq_mod(llm)
+       real hu_mod(llm),hv_mod(llm),w_mod(llm),omega_mod(llm)
+ 
+       integer l,k,k1,k2,kp
+       real aa,frac,frac1,frac2,fact
+ 
+       do l = 1, llm
+
+        if (play(l).ge.plev_prof(nlev_dice)) then
+ 
+        mxcalc=l
+         k1=0
+         k2=0
+
+         if (play(l).le.plev_prof(1)) then
+
+         do k = 1, nlev_dice-1
+          if (play(l).le.plev_prof(k) .and. play(l).gt.plev_prof(k+1)) then
+            k1=k
+            k2=k+1
+          endif
+         enddo
+
+         if (k1.eq.0 .or. k2.eq.0) then
+          write(*,*) 'PB! k1, k2 = ',k1,k2
+          write(*,*) 'l,play(l) = ',l,play(l)/100
+         do k = 1, nlev_dice-1
+          write(*,*) 'k,plev_prof(k) = ',k,plev_prof(k)/100
+         enddo
+         endif
+
+         frac = (plev_prof(k2)-play(l))/(plev_prof(k2)-plev_prof(k1))
+         th_mod(l)= th_prof(k2) - frac*(th_prof(k2)-th_prof(k1))
+         qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
+         u_mod(l)= u_prof(k2) - frac*(u_prof(k2)-u_prof(k1))
+         v_mod(l)= v_prof(k2) - frac*(v_prof(k2)-v_prof(k1))
+         o3_mod(l)= o3_prof(k2) - frac*(o3_prof(k2)-o3_prof(k1))
+         ht_mod(l)= ht_prof(k2) - frac*(ht_prof(k2)-ht_prof(k1))
+         hq_mod(l)= hq_prof(k2) - frac*(hq_prof(k2)-hq_prof(k1))
+         hu_mod(l)= hu_prof(k2) - frac*(hu_prof(k2)-hu_prof(k1))
+         hv_mod(l)= hv_prof(k2) - frac*(hv_prof(k2)-hv_prof(k1))
+         w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
+         omega_mod(l)= omega_prof(k2) - frac*(omega_prof(k2)-omega_prof(k1))
+     
+         else !play>plev_prof(1)
+
+         k1=1
+         k2=2
+         frac1 = (play(l)-plev_prof(k2))/(plev_prof(k1)-plev_prof(k2))
+         frac2 = (play(l)-plev_prof(k1))/(plev_prof(k1)-plev_prof(k2))
+         th_mod(l)= frac1*th_prof(k1) - frac2*th_prof(k2)
+         qv_mod(l)= frac1*qv_prof(k1) - frac2*qv_prof(k2)
+         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
+         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
+         o3_mod(l)= frac1*o3_prof(k1) - frac2*o3_prof(k2)
+         ht_mod(l)= frac1*ht_prof(k1) - frac2*ht_prof(k2)
+         hq_mod(l)= frac1*hq_prof(k1) - frac2*hq_prof(k2)
+         hu_mod(l)= frac1*hu_prof(k1) - frac2*hu_prof(k2)
+         hv_mod(l)= frac1*hv_prof(k1) - frac2*hv_prof(k2)
+         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
+         omega_mod(l)= frac1*omega_prof(k1) - frac2*omega_prof(k2)
+
+         endif ! play.le.plev_prof(1)
+
+        else ! above max altitude of forcing file
+ 
+!jyg
+         fact=20.*(plev_prof(nlev_dice)-play(l))/plev_prof(nlev_dice) !jyg
+         fact = max(fact,0.)                                           !jyg
+         fact = exp(-fact)                                             !jyg
+         th_mod(l)= th_prof(nlev_dice)                                 !jyg
+         qv_mod(l)= qv_prof(nlev_dice)*fact                            !jyg
+         u_mod(l)= u_prof(nlev_dice)*fact                              !jyg
+         v_mod(l)= v_prof(nlev_dice)*fact                              !jyg
+         o3_mod(l)= o3_prof(nlev_dice)*fact                            !jyg
+         ht_mod(l)= ht_prof(nlev_dice)                                 !jyg
+         hq_mod(l)= hq_prof(nlev_dice)*fact                            !jyg
+         hu_mod(l)= hu_prof(nlev_dice)                                 !jyg
+         hv_mod(l)= hv_prof(nlev_dice)                                 !jyg
+         w_mod(l)= 0.                                                  !jyg
+         omega_mod(l)= 0.                                              !jyg
+ 
+        endif ! play
+ 
+       enddo ! l
+
+!       do l = 1,llm
+!       print *,'t_mod(l),q_mod(l),ht_mod(l),hq_mod(l) ',
+!     $        l,t_mod(l),q_mod(l),ht_mod(l),hq_mod(l)
+!       enddo
+ 
+          return
+          end
+
 !======================================================================
         SUBROUTINE interp_astex_time(day,day1,annee_ref                    &
@@ -2652,4 +2775,121 @@
 
 !======================================================================
+        SUBROUTINE interp_dice_time(day,day1,annee_ref                    &
+     &             ,year_ini_dice,day_ini_dice,nt_dice,dt_dice            &
+     &             ,nlev_dice,shf_dice,lhf_dice,lwup_dice,swup_dice       &
+     &             ,tg_dice,ustar_dice,psurf_dice,ug_dice,vg_dice         &
+     &             ,ht_dice,hq_dice,hu_dice,hv_dice,w_dice,omega_dice     &
+     &             ,shf_prof,lhf_prof,lwup_prof,swup_prof,tg_prof         &
+     &             ,ustar_prof,psurf_prof,ug_prof,vg_prof                 &
+     &             ,ht_prof,hq_prof,hu_prof,hv_prof,w_prof,omega_prof)
+        implicit none
+
+!---------------------------------------------------------------------------------------
+! Time interpolation of a 2D field to the timestep corresponding to day
+!
+! day: current julian day (e.g. 717538.2)
+! day1: first day of the simulation
+! nt_dice: total nb of data in the forcing (e.g. 145 for Dice)
+! dt_dice: total time interval (in sec) between 2 forcing data (e.g. 30min. for Dice)
+!---------------------------------------------------------------------------------------
+
+#include "compar1d.h"
+
+! inputs:
+        integer annee_ref
+        integer nt_dice,nlev_dice
+        integer year_ini_dice
+        real day, day1,day_ini_dice,dt_dice
+        real shf_dice(nt_dice),lhf_dice(nt_dice),lwup_dice(nt_dice)
+        real swup_dice(nt_dice),tg_dice(nt_dice),ustar_dice(nt_dice)
+        real psurf_dice(nt_dice),ug_dice(nt_dice),vg_dice(nt_dice)
+        real ht_dice(nlev_dice,nt_dice),hq_dice(nlev_dice,nt_dice)
+        real hu_dice(nlev_dice,nt_dice),hv_dice(nlev_dice,nt_dice)
+        real w_dice(nlev_dice,nt_dice),omega_dice(nlev_dice,nt_dice)
+! outputs:
+        real tg_prof,shf_prof,lhf_prof,lwup_prof,swup_prof
+        real ustar_prof,psurf_prof,ug_prof,vg_prof
+        real ht_prof(nlev_dice),hq_prof(nlev_dice)
+        real hu_prof(nlev_dice),hv_prof(nlev_dice)
+        real w_prof(nlev_dice),omega_prof(nlev_dice)
+! local:
+        integer it_dice1, it_dice2,k
+        real timeit,time_dice1,time_dice2,frac
+
+
+        if (forcing_type.eq.7) then
+! Check that initial day of the simulation consistent with Dice period:
+       print *,'annee_ref=',annee_ref
+       print *,'day1=',day1
+       print *,'day_ini_dice=',day_ini_dice
+       if (annee_ref.ne.1999) then
+        print*,'Pour Dice, annee_ref doit etre 1999'
+        print*,'Changer annee_ref dans run.def'
+        stop
+       endif
+       if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice) then
+        print*,'Dice a debute le 23 Oct 1999 (jour julien=296)'
+        print*,'Changer dayref dans run.def',day1,day_ini_dice
+        stop
+       endif
+       if (annee_ref.eq.1999 .and. day1.gt.day_ini_dice+2) then
+        print*,'Dice a fini le 25 Oct 1999 (jour julien=298)'
+        print*,'Changer dayref ou nday dans run.def',day1,day_ini_dice
+        stop
+       endif
+
+       endif
+
+! Determine timestep relative to the 1st day of TOGA-COARE:
+!       timeit=(day-day1)*86400.
+!       if (annee_ref.eq.1992) then
+!        timeit=(day-day_ini_dice)*86400.
+!       else
+!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
+!       endif
+      timeit=(day-day_ini_dice)*86400
+
+! Determine the closest observation times:
+       it_dice1=INT(timeit/dt_dice)+1
+       it_dice2=it_dice1 + 1
+       time_dice1=(it_dice1-1)*dt_dice
+       time_dice2=(it_dice2-1)*dt_dice
+
+       if (it_dice1 .ge. nt_dice) then
+        write(*,*) 'PB-stop: day, it_dice1, it_dice2, timeit: ',day,it_dice1,it_dice2,timeit/86400.
+        stop
+       endif
+
+! time interpolation:
+       frac=(time_dice2-timeit)/(time_dice2-time_dice1)
+       frac=max(frac,0.0)
+
+       shf_prof = shf_dice(it_dice2)-frac*(shf_dice(it_dice2)-shf_dice(it_dice1))
+       lhf_prof = lhf_dice(it_dice2)-frac*(lhf_dice(it_dice2)-lhf_dice(it_dice1))
+       lwup_prof = lwup_dice(it_dice2)-frac*(lwup_dice(it_dice2)-lwup_dice(it_dice1))
+       swup_prof = swup_dice(it_dice2)-frac*(swup_dice(it_dice2)-swup_dice(it_dice1))
+       tg_prof = tg_dice(it_dice2)-frac*(tg_dice(it_dice2)-tg_dice(it_dice1))
+       ustar_prof = ustar_dice(it_dice2)-frac*(ustar_dice(it_dice2)-ustar_dice(it_dice1))
+       psurf_prof = psurf_dice(it_dice2)-frac*(psurf_dice(it_dice2)-psurf_dice(it_dice1))
+       ug_prof = ug_dice(it_dice2)-frac*(ug_dice(it_dice2)-ug_dice(it_dice1))
+       vg_prof = vg_dice(it_dice2)-frac*(vg_dice(it_dice2)-vg_dice(it_dice1))
+
+!        print*,
+!     :'day,annee_ref,day_ini_dice,timeit,it_dice1,it_dice2,SST:',
+!     :day,annee_ref,day_ini_dice,timeit/86400.,it_dice1,it_dice2,ts_prof
+
+       do k=1,nlev_dice
+        ht_prof(k) = ht_dice(k,it_dice2)-frac*(ht_dice(k,it_dice2)-ht_dice(k,it_dice1))
+        hq_prof(k) = hq_dice(k,it_dice2)-frac*(hq_dice(k,it_dice2)-hq_dice(k,it_dice1))
+        hu_prof(k) = hu_dice(k,it_dice2)-frac*(hu_dice(k,it_dice2)-hu_dice(k,it_dice1))
+        hv_prof(k) = hv_dice(k,it_dice2)-frac*(hv_dice(k,it_dice2)-hv_dice(k,it_dice1))
+        w_prof(k) = w_dice(k,it_dice2)-frac*(w_dice(k,it_dice2)-w_dice(k,it_dice1))
+        omega_prof(k) = omega_dice(k,it_dice2)-frac*(omega_dice(k,it_dice2)-omega_dice(k,it_dice1))
+        enddo
+
+        return
+        END
+
+!======================================================================
         SUBROUTINE interp_armcu_time(day,day1,annee_ref                    &
      &             ,year_ini_armcu,day_ini_armcu,nt_armcu,dt_armcu         &
@@ -2916,341 +3156,4 @@
         return
         end
-!=====================================================================
-      subroutine read_amma(fich_amma,nlevel,ntime                          &
-     &     ,zz,pp,temp,qv,u,v,dw                                           &
-     &     ,dt,dq,sens,flat)
-
-!program reading forcings of the AMMA case study
-
-
-      implicit none
-
-#include "netcdf.inc"
-
-      integer ntime,nlevel
-      character*80 :: fich_amma
-      real*8 zz(nlevel)
-
-      real*8 temp(nlevel),pp(nlevel)
-      real*8 qv(nlevel),u(nlevel)
-      real*8 v(nlevel)
-      real*8 dw(nlevel,ntime)
-      real*8 dt(nlevel,ntime)
-      real*8 dq(nlevel,ntime)
-      real*8 flat(ntime),sens(ntime)
-
-      integer nid, ierr
-      integer nbvar3d
-      parameter(nbvar3d=30)
-      integer var3didin(nbvar3d)
-
-      ierr = NF_OPEN(fich_amma,NF_NOWRITE,nid)
-      if (ierr.NE.NF_NOERR) then
-         write(*,*) 'ERROR: Pb opening forcings nc file '
-         write(*,*) NF_STRERROR(ierr)
-         stop ""
-      endif
-
-
-       ierr=NF_INQ_VARID(nid,"zz",var3didin(1)) 
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'lev'
-         endif
-
-
-      ierr=NF_INQ_VARID(nid,"temp",var3didin(2))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'temp'
-         endif
-
-      ierr=NF_INQ_VARID(nid,"qv",var3didin(3))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'qv'
-         endif
-
-      ierr=NF_INQ_VARID(nid,"u",var3didin(4))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'u'
-         endif
-
-      ierr=NF_INQ_VARID(nid,"v",var3didin(5))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'v'
-         endif
-
-      ierr=NF_INQ_VARID(nid,"dw",var3didin(6))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'dw'
-         endif
-
-      ierr=NF_INQ_VARID(nid,"dt",var3didin(7))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'dt'
-         endif
-
-      ierr=NF_INQ_VARID(nid,"dq",var3didin(8))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'dq'
-         endif
-      
-      ierr=NF_INQ_VARID(nid,"sens",var3didin(9))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'sens'
-         endif
-
-      ierr=NF_INQ_VARID(nid,"flat",var3didin(10))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'flat'
-         endif
-
-      ierr=NF_INQ_VARID(nid,"pp",var3didin(11))
-         if(ierr/=NF_NOERR) then
-           write(*,*) NF_STRERROR(ierr)
-           stop 'pp'
-      endif
-
-!dimensions lecture
-!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
- 
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture z ok',zz
-
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),temp)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(2),temp)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture th ok',temp
-
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qv)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(3),qv)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture qv ok',qv
- 
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),u)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(4),u)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture u ok',u
-
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),v)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(5),v)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture v ok',v
-
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),dw)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(6),dw)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture w ok',dw
-
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),dt)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(7),dt)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture dt ok',dt
-
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),dq)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(8),dq)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture dq ok',dq
-
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),sens)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(9),sens)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture sens ok',sens
-
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),flat)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(10),flat)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture flat ok',flat
-
-#ifdef NC_DOUBLE
-         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pp)
-#else
-         ierr = NF_GET_VAR_REAL(nid,var3didin(11),pp)
-#endif
-         if(ierr/=NF_NOERR) then
-            write(*,*) NF_STRERROR(ierr)
-            stop "getvarup"
-         endif
-!          write(*,*)'lecture pp ok',pp
-
-         return 
-         end subroutine read_amma
-!======================================================================
-        SUBROUTINE interp_amma_time(day,day1,annee_ref                     &
-     &         ,year_ini_amma,day_ini_amma,nt_amma,dt_amma,nlev_amma       &
-     &         ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma               &
-     &         ,vitw_prof,ht_prof,hq_prof,lat_prof,sens_prof)
-        implicit none
-
-!---------------------------------------------------------------------------------------
-! Time interpolation of a 2D field to the timestep corresponding to day
-!
-! day: current julian day (e.g. 717538.2)
-! day1: first day of the simulation
-! nt_amma: total nb of data in the forcing (e.g. 48 for AMMA)
-! dt_amma: total time interval (in sec) between 2 forcing data (e.g. 30min for AMMA)
-!---------------------------------------------------------------------------------------
-
-#include "compar1d.h"
-
-! inputs:
-        integer annee_ref
-        integer nt_amma,nlev_amma
-        integer year_ini_amma
-        real day, day1,day_ini_amma,dt_amma
-        real vitw_amma(nlev_amma,nt_amma)
-        real ht_amma(nlev_amma,nt_amma)
-        real hq_amma(nlev_amma,nt_amma)
-        real lat_amma(nt_amma)
-        real sens_amma(nt_amma)
-! outputs:
-        real vitw_prof(nlev_amma)
-        real ht_prof(nlev_amma)
-        real hq_prof(nlev_amma)
-        real lat_prof,sens_prof
-! local:
-        integer it_amma1, it_amma2,k
-        real timeit,time_amma1,time_amma2,frac
-
-
-        if (forcing_type.eq.6) then
-! Check that initial day of the simulation consistent with AMMA case:
-       if (annee_ref.ne.2006) then
-        print*,'Pour AMMA, annee_ref doit etre 2006'
-        print*,'Changer annee_ref dans run.def'
-        stop
-       endif
-       if (annee_ref.eq.2006 .and. day1.lt.day_ini_amma) then
-        print*,'AMMA a débuté le 10 juillet 2006',day1,day_ini_amma
-        print*,'Changer dayref dans run.def'
-        stop
-       endif
-       if (annee_ref.eq.2006 .and. day1.gt.day_ini_amma+1) then
-        print*,'AMMA a fini le 11 juillet'
-        print*,'Changer dayref ou nday dans run.def'
-        stop
-       endif
-       endif
-
-! Determine timestep relative to the 1st day of AMMA:
-!       timeit=(day-day1)*86400.
-!       if (annee_ref.eq.1992) then
-!        timeit=(day-day_ini_toga)*86400.
-!       else
-!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
-!       endif
-      timeit=(day-day_ini_amma)*86400
-
-! Determine the closest observation times:
-!       it_amma1=INT(timeit/dt_amma)+1
-!       it_amma2=it_amma1 + 1
-!       time_amma1=(it_amma1-1)*dt_amma
-!       time_amma2=(it_amma2-1)*dt_amma
-
-       it_amma1=INT(timeit/dt_amma)+1
-       IF (it_amma1 .EQ. nt_amma) THEN
-       it_amma2=it_amma1 
-       ELSE
-       it_amma2=it_amma1 + 1
-       ENDIF
-       time_amma1=(it_amma1-1)*dt_amma
-       time_amma2=(it_amma2-1)*dt_amma
-
-       if (it_amma1 .gt. nt_amma) then
-        write(*,*) 'PB-stop: day, it_amma1, it_amma2, timeit: '            &
-     &        ,day,day_ini_amma,it_amma1,it_amma2,timeit/86400.
-        stop
-       endif
-
-! time interpolation:
-       frac=(time_amma2-timeit)/(time_amma2-time_amma1)
-       frac=max(frac,0.0)
-
-       lat_prof = lat_amma(it_amma2)                                       &
-     &          -frac*(lat_amma(it_amma2)-lat_amma(it_amma1)) 
-       sens_prof = sens_amma(it_amma2)                                     &
-     &          -frac*(sens_amma(it_amma2)-sens_amma(it_amma1))
-
-       do k=1,nlev_amma
-        vitw_prof(k) = vitw_amma(k,it_amma2)                               &
-     &          -frac*(vitw_amma(k,it_amma2)-vitw_amma(k,it_amma1))
-        ht_prof(k) = ht_amma(k,it_amma2)                                   &
-     &          -frac*(ht_amma(k,it_amma2)-ht_amma(k,it_amma1))
-        hq_prof(k) = hq_amma(k,it_amma2)                                   &
-     &          -frac*(hq_amma(k,it_amma2)-hq_amma(k,it_amma1))
-        enddo
-
-        return
-        END
 
 !=====================================================================
@@ -3516,2 +3419,563 @@
          return 
          end subroutine read_fire
+!=====================================================================
+      subroutine read_dice(fich_dice,nlevel,ntime                         &
+     &     ,zz,pres,th,qv,u,v,o3                                          &
+     &     ,shf,lhf,lwup,swup,tg,ustar,psurf,ug,vg                        &
+     &     ,hadvt,hadvq,hadvu,hadvv,w,omega)
+
+!program reading initial profils and forcings of the Dice case study
+
+
+      implicit none
+
+#include "netcdf.inc"
+
+      integer ntime,nlevel
+      integer l,k
+      character*80 :: fich_dice
+      real*8 time(ntime)
+      real*8 zz(nlevel)
+
+      real*8 th(nlevel),pres(nlevel)
+      real*8 qv(nlevel),u(nlevel),v(nlevel),o3(nlevel)
+      real*8 shf(ntime),lhf(ntime),lwup(ntime),swup(ntime),tg(ntime)
+      real*8 ustar(ntime),psurf(ntime),ug(ntime),vg(ntime)
+      real*8 hadvt(nlevel,ntime),hadvq(nlevel,ntime),hadvu(nlevel,ntime)
+      real*8 hadvv(nlevel,ntime),w(nlevel,ntime),omega(nlevel,ntime)
+
+      integer nid, ierr
+      integer nbvar3d
+      parameter(nbvar3d=30)
+      integer var3didin(nbvar3d)
+
+      ierr = NF_OPEN(fich_dice,NF_NOWRITE,nid)
+      if (ierr.NE.NF_NOERR) then
+         write(*,*) 'ERROR: Pb opening forcings nc file '
+         write(*,*) NF_STRERROR(ierr)
+         stop ""
+      endif
+
+
+       ierr=NF_INQ_VARID(nid,"height",var3didin(1)) 
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'height'
+         endif
+
+       ierr=NF_INQ_VARID(nid,"pf",var3didin(11)) 
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'pf'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"theta",var3didin(12))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'theta'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"qv",var3didin(13))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'qv'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"u",var3didin(14))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'u'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"v",var3didin(15))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'v'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"o3mmr",var3didin(16))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'o3'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"shf",var3didin(2))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'shf'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"lhf",var3didin(3))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'lhf'
+         endif
+      
+      ierr=NF_INQ_VARID(nid,"lwup",var3didin(4))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'lwup'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"swup",var3didin(5))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'dqtdx'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"Tg",var3didin(6))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'Tg'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"ustar",var3didin(7))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'ustar'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"psurf",var3didin(8))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'psurf'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"Ug",var3didin(9))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'Ug'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"Vg",var3didin(10))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'Vg'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"hadvT",var3didin(17))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'hadvT'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"hadvq",var3didin(18))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'hadvq'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"hadvu",var3didin(19))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'hadvu'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"hadvv",var3didin(20))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'hadvv'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"w",var3didin(21))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'w'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"omega",var3didin(22))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'omega'
+      endif
+!dimensions lecture
+!      call catchaxis(nid,ntime,nlevel,time,z,ierr)
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(1),zz)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(1),zz)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture zz ok',zz
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),pres)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(11),pres)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture pres ok',pres
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),th)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(12),th)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture th ok',th
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),qv)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(13),qv)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture qv ok',qv
+ 
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(14),u)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(14),u)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture u ok',u
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(15),v)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(15),v)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture v ok',v
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(16),o3)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(16),o3)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture o3 ok',o3
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(2),shf)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(2),shf)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture shf ok',shf
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),lhf)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(3),lhf)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture lhf ok',lhf
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(4),lwup)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(4),lwup)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture lwup ok',lwup
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(5),swup)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(5),swup)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture swup ok',swup
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(6),tg)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(6),tg)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture tg ok',tg
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ustar)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(7),ustar)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture ustar ok',ustar
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(8),psurf)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(8),psurf)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture psurf ok',psurf
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(9),ug)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(9),ug)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture ug ok',ug
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),vg)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(10),vg)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture vg ok',vg
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(17),hadvt)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(17),hadvt)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture hadvt ok',hadvt
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(18),hadvq)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(18),hadvq)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture hadvq ok',hadvq
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(19),hadvu)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(19),hadvu)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture hadvu ok',hadvu
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(20),hadvv)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(20),hadvv)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture hadvv ok',hadvv
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(21),w)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(21),w)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture w ok',w
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(22),omega)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(22),omega)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture omega ok',omega
+
+         return 
+         end subroutine read_dice
+!=====================================================================
+!     Reads CIRC input files      
+
+      SUBROUTINE read_circ(nlev_circ,cf,lwp,iwp,reliq,reice,t,z,p,pm,h2o,o3,sza)
+      
+      parameter (ncm_1=49180)
+#include "YOMCST.h"
+
+      real albsfc(ncm_1), albsfc_w(ncm_1)
+      real cf(nlev_circ), icefra(nlev_circ), deice(nlev_circ), &
+           reliq(nlev_circ), reice(nlev_circ), lwp(nlev_circ), iwp(nlev_circ)
+      real t(nlev_circ+1), z(nlev_circ+1), dz(nlev_circ), p(nlev_circ+1)
+      real aer_beta(nlev_circ), waer(nlev_circ), gaer(nlev_circ)
+      real pm(nlev_circ), tm(nlev_circ), h2o(nlev_circ), o3(nlev_circ)
+      real co2(nlev_circ), n2o(nlev_circ), co(nlev_circ), ch4(nlev_circ), &
+           o2(nlev_circ), ccl4(nlev_circ), f11(nlev_circ), f12(nlev_circ)
+!     za= zenital angle
+!     sza= cosinus angle zenital
+      real wavn(ncm_1), ssf(ncm_1),za,sza
+      integer nlev
+
+
+!     Open the files
+
+      open (11, file='Tsfc_sza_nlev_case.txt', status='old')
+      open (12, file='level_input_case.txt', status='old')
+      open (13, file='layer_input_case.txt', status='old')
+      open (14, file='aerosol_input_case.txt', status='old')
+      open (15, file='cloud_input_case.txt', status='old')
+      open (16, file='sfcalbedo_input_case.txt', status='old')
+      
+!     Read scalar information
+      do iskip=1,5
+         read (11, *)
+      enddo
+      read (11, '(i8)') nlev
+      read (11, '(f10.2)') tsfc
+      read (11, '(f10.2)') za
+      read (11, '(f10.4)') sw_dn_toa
+      sza=cos(za/180.*RPI)
+      print *,'nlev,tsfc,sza,sw_dn_toa,RPI',nlev,tsfc,sza,sw_dn_toa,RPI
+      close(11)
+
+!     Read level information
+      read (12, *)
+      do il=1,nlev
+         read (12, 302) ilev, z(il), p(il), t(il)
+         z(il)=z(il)*1000.    ! z donne en km
+         p(il)=p(il)*100.     ! p donne en mb
+      enddo
+302   format (i8, f8.3, 2f9.2)
+      close(12)
+
+!     Read layer information (midpoint values)
+      do iskip=1,3
+         read (13, *)
+      enddo
+      do il=1,nlev-1
+         read (13, 303) ilev,pm(il),tm(il),h2o(il),co2(il),o3(il), &
+                        n2o(il),co(il),ch4(il),o2(il),ccl4(il), &
+                        f11(il),f12(il)
+         pm(il)=pm(il)*100.
+      enddo
+303   format (i8, 2f9.2, 10(2x,e13.7))      
+      close(13)
+      
+!     Read aerosol layer information
+      do iskip=1,3
+         read (14, *)
+      enddo
+      read (14, '(f10.2)') aer_alpha
+      read (14, *)
+      read (14, *)
+      do il=1,nlev-1
+         read (14, 304) ilev, aer_beta(il), waer(il), gaer(il)
+      enddo
+304   format (i8, f9.5, 2f8.3)
+      close(14)
+      
+!     Read cloud information
+      do iskip=1,3
+         read (15, *)
+      enddo
+      do il=1,nlev-1
+         read (15, 305) ilev, cf(il), lwp(il), iwp(il), reliq(il), reice(il)
+         lwp(il)=lwp(il)/1000.          ! lwp donne en g/kg
+         iwp(il)=iwp(il)/1000.          ! iwp donne en g/kg
+         reliq(il)=reliq(il)/1000000.   ! reliq donne en microns
+         reice(il)=reice(il)/1000000.   ! reice donne en microns
+      enddo
+305   format (i8, f8.3, 4f9.2)
+      close(15)
+
+!     Read surface albedo (weighted & unweighted) and spectral solar irradiance
+      do iskip=1,6
+         read (16, *)
+      enddo
+      do icm_1=1,ncm_1
+         read (16, 306) wavn(icm_1), albsfc(icm_1), albsfc_w(icm_1), ssf(icm_1)
+      enddo
+306   format(f10.1, 2f12.5, f14.8)
+      close(16)
+ 
+      return 
+      end subroutine read_circ
+!=====================================================================
+!     Reads RTMIP input files      
+
+      SUBROUTINE read_rtmip(nlev_rtmip,play,plev,t,h2o,o3)
+      
+#include "YOMCST.h"
+
+      real t(nlev_rtmip), pt(nlev_rtmip),pb(nlev_rtmip),h2o(nlev_rtmip), o3(nlev_rtmip)
+      real temp(nlev_rtmip), play(nlev_rtmip),ovap(nlev_rtmip), oz(nlev_rtmip),plev(nlev_rtmip+1)
+      integer nlev
+
+
+!     Open the files
+
+      open (11, file='low_resolution_profile.txt', status='old')
+      
+!     Read level information
+      read (11, *)
+      do il=1,nlev_rtmip
+         read (11, 302) pt(il), pb(il), t(il),h2o(il),o3(il)
+      enddo
+      do il=1,nlev_rtmip
+         play(il)=pt(nlev_rtmip-il+1)*100.     ! p donne en mb
+         temp(il)=t(nlev_rtmip-il+1)
+         ovap(il)=h2o(nlev_rtmip-il+1)
+         oz(il)=o3(nlev_rtmip-il+1)
+      enddo
+      do il=1,39
+         plev(il)=play(il)+(play(il+1)-play(il))/2.
+         print *,'il p t ovap oz=',il,plev(il),temp(il),ovap(il),oz(il)
+      enddo
+      plev(41)=101300.
+302   format (e16.10,3x,e16.10,3x,e16.10,3x,e12.6,3x,e12.6)
+      close(12)
+ 
+      return 
+      end subroutine read_rtmip
+!=====================================================================
+
