Changeset 1795 for LMDZ5/branches/testing/libf/phy1d/1DUTILS.h_no_writelim

Timestamp:

Jul 18, 2013, 10:20:28 AM (11 years ago)

Author:

Ehouarn Millour

Message:

Version testing basee sur la r1794

---

Testing release based on r1794

Location:

LMDZ5/branches/testing

Files:

• . (modified) (1 prop)
• libf/phy1d/1DUTILS.h_no_writelim (modified) (13 diffs)

LMDZ5/branches/testing/libf/phy1d/1DUTILS.h_no_writelim

@@ -27 +27 @@
 #include "compar1d.h"
 #include "flux_arp.h"
+#include "tsoilnudge.h"
 #include "fcg_gcssold.h"
 #include "fcg_racmo.h"
@@ -100 +101 @@
 !             initial profiles from RICO idealized
 !             LS convergence imposed from  RICO (cst) 
+!         = 6 ==> forcing_amma = .true.
 !         = 40 ==> forcing_GCSSold = .true.
 !             initial profile from GCSS file
 !             LS convergence imposed from GCSS file
+!         = 59 ==> forcing_sandu = .true.
+!             initial profiles from sanduref file: see prof.inp.001
+!             SST varying with time and divergence constante: see ifa_sanduref.txt file
+!             Radiation has to be computed interactively
+!         = 60 ==> forcing_astex = .true.
+!             initial profiles from file: see prof.inp.001 
+!             SST,divergence,ug,vg,ufa,vfa varying with time : see ifa_astex.txt file
+!             Radiation has to be computed interactively
+!         = 61 ==> forcing_armcu = .true.
+!             initial profiles from file: see prof.inp.001 
+!             sensible and latent heat flux imposed: see ifa_arm_cu_1.txt
+!             large scale advective forcing & radiative tendencies applied below 1000m: see ifa_arm_cu_2.txt
+!             use geostrophic wind ug=10m/s vg=0m/s. Duration of the case 53100s 
+!             Radiation to be switched off
 !
        forcing_type = 0
@@ -126 +142 @@
        CALL getin('ok_flux_surf',ok_flux_surf)
 
+!Config  Key  = ok_old_disvert
+!Config  Desc = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
+!Config  Def  = false
+!Config  Help = utilisation de l ancien programme disvert0 (dans 1DUTILS.h)
+       ok_old_disvert = .FALSE.
+       CALL getin('ok_old_disvert',ok_old_disvert)
+
 !Config  Key  = time_ini
 !Config  Desc = meaningless in this  case
@@ -227 +250 @@
        zpicinp = 300.
        CALL getin('zpicinp',zpicinp)
+!Config key = nudge_tsoil
+!Config  Desc = activation of soil temperature nudging
+!Config  Def  = .FALSE.
+!Config  Help = ...
+
+       nudge_tsoil=.FALSE.
+       CALL getin('nudge_tsoil',nudge_tsoil)
+
+!Config key = isoil_nudge
+!Config  Desc = level number where soil temperature is nudged
+!Config  Def  = 3
+!Config  Help = ...
+
+       isoil_nudge=3
+       CALL getin('isoil_nudge',isoil_nudge)
+
+!Config key = Tsoil_nudge
+!Config  Desc = target temperature for tsoil(isoil_nudge)
+!Config  Def  = 300.
+!Config  Help = ...
+
+       Tsoil_nudge=300.
+       CALL getin('Tsoil_nudge',Tsoil_nudge)
+
+!Config key = tau_soil_nudge
+!Config  Desc = nudging relaxation time for tsoil
+!Config  Def  = 3600.
+!Config  Help = ...
+
+       tau_soil_nudge=3600.
+       CALL getin('tau_soil_nudge',tau_soil_nudge)
+
+
 
 
@@ -250 +306 @@
       write(lunout,*)' qsolinp = ', qsolinp
       write(lunout,*)' zpicinp = ', zpicinp
+      write(lunout,*)' nudge_tsoil = ', nudge_tsoil
+      write(lunout,*)' isoil_nudge = ', isoil_nudge
+      write(lunout,*)' Tsoil_nudge = ', Tsoil_nudge
+      write(lunout,*)' tau_soil_nudge = ', tau_soil_nudge
       IF (forcing_type .eq.40) THEN
         write(lunout,*) '--- Forcing type GCSS Old --- with:'
@@ -703 +763 @@
       RETURN
       END
-      subroutine scopy(n,sx,incx,sy,incy)
-!
-      IMPLICIT NONE
-!
-      integer n,incx,incy,ix,iy,i
-      real sx((n-1)*incx+1),sy((n-1)*incy+1)
-!
-      iy=1
-      ix=1
-      do 10 i=1,n
-         sy(iy)=sx(ix)
-         ix=ix+incx
-         iy=iy+incy
-10    continue
-!
-      return
-      end
       subroutine wrgradsfi(if,nl,field,name,titlevar)
       implicit none
@@ -956 +999 @@
       END
  
-      SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
- 
+      SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
+ 
+!    Ancienne version disvert dont on a modifie nom pour utiliser
+!    le disvert de dyn3d (qui permet d'utiliser grille avec ab,bp imposes)
+!    (MPL 18092012)
+!
 !    Auteur :  P. Le Van .
 !
@@ -1402 +1449 @@
           end
 
+c-------------------------------------------------------------------------
+      SUBROUTINE read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
+      implicit none
+
+c-------------------------------------------------------------------------
+c Read I.SANDU case forcing data
+c-------------------------------------------------------------------------
+
+      integer nlev_sandu,nt_sandu
+      real ts_sandu(nt_sandu)
+      character*80 fich_sandu
+
+      integer no,l,k,ip
+      real riy,rim,rid,rih,bid
+
+      integer iy,im,id,ih
+
+       real plev_min
+
+       plev_min = 55000.  ! pas de tendance de vap. d eau au-dessus de 55 hPa
+
+      open(21,file=trim(fich_sandu),form='formatted')
+      read(21,'(a)')
+      do ip = 1, nt_sandu
+      read(21,'(a)')
+      read(21,'(a)')
+      read(21,223) iy, im, id, ih, ts_sandu(ip)
+      print *,'ts=',iy,im,id,ih,ip,ts_sandu(ip)
+      enddo
+      close(21)
+
+  223 format(4i3,f8.2)
+  226 format(f7.1,1x,10f8.2)
+  227 format(f7.1,1x,1p,4e11.3)
+  230 format(6f9.3,4e11.3)
+
+          return
+          end
+
+!=====================================================================
+c-------------------------------------------------------------------------
+      SUBROUTINE read_astex(fich_astex,nlev_astex,nt_astex,div_astex,
+     : ts_astex,ug_astex,vg_astex,ufa_astex,vfa_astex)
+      implicit none
+
+c-------------------------------------------------------------------------
+c Read Astex case forcing data
+c-------------------------------------------------------------------------
+
+      integer nlev_astex,nt_astex
+      real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
+      real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
+      character*80 fich_astex
+
+      integer no,l,k,ip
+      real riy,rim,rid,rih,bid
+
+      integer iy,im,id,ih
+
+       real plev_min
+
+       plev_min = 55000.  ! pas de tendance de vap. d eau au-dessus de 55 hPa
+
+      open(21,file=trim(fich_astex),form='formatted')
+      read(21,'(a)')
+      read(21,'(a)')
+      do ip = 1, nt_astex
+      read(21,'(a)')
+      read(21,'(a)')
+      read(21,223) iy, im, id, ih, div_astex(ip),ts_astex(ip),
+     :ug_astex(ip),vg_astex(ip),ufa_astex(ip),vfa_astex(ip)
+      ts_astex(ip)=ts_astex(ip)+273.15
+      print *,'ts=',iy,im,id,ih,ip,div_astex(ip),ts_astex(ip),
+     :ug_astex(ip),vg_astex(ip),ufa_astex(ip),vg_astex(ip)
+      enddo
+      close(21)
+
+  223 format(4i3,e13.2,f7.2,f7.3,f7.2,f7.3,f7.2)
+  226 format(f7.1,1x,10f8.2)
+  227 format(f7.1,1x,1p,4e11.3)
+  230 format(6f9.3,4e11.3)
+
+          return
+          end
 !=====================================================================
       subroutine read_twpice(fich_twpice,nlevel,ntime
@@ -1884 +2015 @@
        return
        end
+!=====================================================================
+
+       SUBROUTINE interp_sandu_vertical(play,nlev_sandu,plev_prof
+     :         ,t_prof,thl_prof,q_prof,u_prof,v_prof,w_prof
+     :         ,omega_prof,o3mmr_prof
+     :         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod
+     :         ,omega_mod,o3mmr_mod,mxcalc)
+
+       implicit none
+
+#include "dimensions.h"
+
+c-------------------------------------------------------------------------
+c Vertical interpolation of SANDUREF forcing data onto model levels
+c-------------------------------------------------------------------------
+
+       integer nlevmax
+       parameter (nlevmax=41)
+       integer nlev_sandu,mxcalc
+!       real play(llm), plev_prof(nlevmax)
+!       real t_prof(nlevmax),q_prof(nlevmax)
+!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
+!       real ht_prof(nlevmax),vt_prof(nlevmax)
+!       real hq_prof(nlevmax),vq_prof(nlevmax)
+
+       real play(llm), plev_prof(nlev_sandu)
+       real t_prof(nlev_sandu),thl_prof(nlev_sandu),q_prof(nlev_sandu)
+       real u_prof(nlev_sandu),v_prof(nlev_sandu), w_prof(nlev_sandu)
+       real omega_prof(nlev_sandu),o3mmr_prof(nlev_sandu)
+
+       real t_mod(llm),thl_mod(llm),q_mod(llm)
+       real u_mod(llm),v_mod(llm), w_mod(llm)
+       real omega_mod(llm),o3mmr_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_sandu)) then
+
+        mxcalc=l
+         k1=0
+         k2=0
+
+         if (play(l).le.plev_prof(1)) then
+
+         do k = 1, nlev_sandu-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_sandu-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))
+         t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
+         thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
+         q_mod(l)= q_prof(k2) - frac*(q_prof(k2)-q_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))
+         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))
+         o3mmr_mod(l)=o3mmr_prof(k2)-
+     :      frac*(o3mmr_prof(k2)-o3mmr_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))
+         t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)
+         thl_mod(l)= frac1*thl_prof(k1) - frac2*thl_prof(k2)
+         q_mod(l)= frac1*q_prof(k1) - frac2*q_prof(k2)
+         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
+         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
+         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
+         omega_mod(l)= frac1*omega_prof(k1) - frac2*omega_prof(k2)
+         o3mmr_mod(l)= frac1*o3mmr_prof(k1) - frac2*o3mmr_prof(k2)
+
+         endif ! play.le.plev_prof(1)
+
+        else ! above max altitude of forcing file
+
+cjyg
+         fact=20.*(plev_prof(nlev_sandu)-play(l))/plev_prof(nlev_sandu) !jyg
+         fact = max(fact,0.)                                           !jyg
+         fact = exp(-fact)                                             !jyg
+         t_mod(l)= t_prof(nlev_sandu)                                   !jyg
+         thl_mod(l)= thl_prof(nlev_sandu)                                   !jyg
+         q_mod(l)= q_prof(nlev_sandu)*fact                              !jyg
+         u_mod(l)= u_prof(nlev_sandu)*fact                              !jyg
+         v_mod(l)= v_prof(nlev_sandu)*fact                              !jyg
+         w_mod(l)= w_prof(nlev_sandu)*fact                              !jyg
+         omega_mod(l)= omega_prof(nlev_sandu)*fact                      !jyg
+         o3mmr_mod(l)= o3mmr_prof(nlev_sandu)*fact                      !jyg
+
+        endif ! play
+
+       enddo ! l
+
+       do l = 1,llm
+!      print *,'t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l) ',
+!    $        l,t_mod(l),thl_mod(l),q_mod(l),u_mod(l),v_mod(l)
+       enddo
+
+          return
+          end
+!=====================================================================
+       SUBROUTINE interp_astex_vertical(play,nlev_astex,plev_prof
+     :         ,t_prof,thl_prof,qv_prof,ql_prof,qt_prof,u_prof,v_prof
+     :         ,w_prof,tke_prof,o3mmr_prof
+     :         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod
+     :         ,tke_mod,o3mmr_mod,mxcalc)
+
+       implicit none
+
+#include "dimensions.h"
+
+c-------------------------------------------------------------------------
+c Vertical interpolation of Astex forcing data onto model levels
+c-------------------------------------------------------------------------
+
+       integer nlevmax
+       parameter (nlevmax=41)
+       integer nlev_astex,mxcalc
+!       real play(llm), plev_prof(nlevmax)
+!       real t_prof(nlevmax),qv_prof(nlevmax)
+!       real u_prof(nlevmax),v_prof(nlevmax), w_prof(nlevmax)
+!       real ht_prof(nlevmax),vt_prof(nlevmax)
+!       real hq_prof(nlevmax),vq_prof(nlevmax)
+
+       real play(llm), plev_prof(nlev_astex)
+       real t_prof(nlev_astex),thl_prof(nlev_astex),qv_prof(nlev_astex)
+       real u_prof(nlev_astex),v_prof(nlev_astex), w_prof(nlev_astex)
+       real o3mmr_prof(nlev_astex),ql_prof(nlev_astex)
+       real qt_prof(nlev_astex),tke_prof(nlev_astex)
+
+       real t_mod(llm),thl_mod(llm),qv_mod(llm)
+       real u_mod(llm),v_mod(llm), w_mod(llm),tke_mod(llm)
+       real o3mmr_mod(llm),ql_mod(llm),qt_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_astex)) then
+
+        mxcalc=l
+         k1=0
+         k2=0
+
+         if (play(l).le.plev_prof(1)) then
+
+         do k = 1, nlev_astex-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_astex-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))
+         t_mod(l)= t_prof(k2) - frac*(t_prof(k2)-t_prof(k1))
+         thl_mod(l)= thl_prof(k2) - frac*(thl_prof(k2)-thl_prof(k1))
+         qv_mod(l)= qv_prof(k2) - frac*(qv_prof(k2)-qv_prof(k1))
+         ql_mod(l)= ql_prof(k2) - frac*(ql_prof(k2)-ql_prof(k1))
+         qt_mod(l)= qt_prof(k2) - frac*(qt_prof(k2)-qt_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))
+         w_mod(l)= w_prof(k2) - frac*(w_prof(k2)-w_prof(k1))
+         tke_mod(l)= tke_prof(k2) - frac*(tke_prof(k2)-tke_prof(k1))
+         o3mmr_mod(l)=o3mmr_prof(k2)-
+     :      frac*(o3mmr_prof(k2)-o3mmr_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))
+         t_mod(l)= frac1*t_prof(k1) - frac2*t_prof(k2)
+         thl_mod(l)= frac1*thl_prof(k1) - frac2*thl_prof(k2)
+         qv_mod(l)= frac1*qv_prof(k1) - frac2*qv_prof(k2)
+         ql_mod(l)= frac1*ql_prof(k1) - frac2*ql_prof(k2)
+         qt_mod(l)= frac1*qt_prof(k1) - frac2*qt_prof(k2)
+         u_mod(l)= frac1*u_prof(k1) - frac2*u_prof(k2)
+         v_mod(l)= frac1*v_prof(k1) - frac2*v_prof(k2)
+         w_mod(l)= frac1*w_prof(k1) - frac2*w_prof(k2)
+         tke_mod(l)= frac1*tke_prof(k1) - frac2*tke_prof(k2)
+         o3mmr_mod(l)= frac1*o3mmr_prof(k1) - frac2*o3mmr_prof(k2)
+
+         endif ! play.le.plev_prof(1)
+
+        else ! above max altitude of forcing file
+
+cjyg
+         fact=20.*(plev_prof(nlev_astex)-play(l))/plev_prof(nlev_astex) !jyg
+         fact = max(fact,0.)                                           !jyg
+         fact = exp(-fact)                                             !jyg
+         t_mod(l)= t_prof(nlev_astex)                                   !jyg
+         thl_mod(l)= thl_prof(nlev_astex)                                   !jyg
+         qv_mod(l)= qv_prof(nlev_astex)*fact                              !jyg
+         ql_mod(l)= ql_prof(nlev_astex)*fact                              !jyg
+         qt_mod(l)= qt_prof(nlev_astex)*fact                              !jyg
+         u_mod(l)= u_prof(nlev_astex)*fact                              !jyg
+         v_mod(l)= v_prof(nlev_astex)*fact                              !jyg
+         w_mod(l)= w_prof(nlev_astex)*fact                              !jyg
+         tke_mod(l)= tke_prof(nlev_astex)*fact                              !jyg
+         o3mmr_mod(l)= o3mmr_prof(nlev_astex)*fact                      !jyg
+
+        endif ! play
+
+       enddo ! l
+
+       do l = 1,llm
+!      print *,'t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l) ',
+!    $        l,t_mod(l),thl_mod(l),qv_mod(l),u_mod(l),v_mod(l)
+       enddo
+
+          return
+          end
 
 !======================================================================
@@ -2048 +2420 @@
           end
 
+!======================================================================
+        SUBROUTINE interp_sandu_time(day,day1,annee_ref
+     i             ,year_ini_sandu,day_ini_sandu,nt_sandu,dt_sandu
+     i             ,nlev_sandu,ts_sandu,ts_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_sandu: total nb of data in the forcing (e.g. 13 for Sanduref)
+! dt_sandu: total time interval (in sec) between 2 forcing data (e.g. 6h for Sanduref)
+!---------------------------------------------------------------------------------------
+! inputs:
+        integer annee_ref
+        integer nt_sandu,nlev_sandu
+        integer year_ini_sandu
+        real day, day1,day_ini_sandu,dt_sandu
+        real ts_sandu(nt_sandu)
+! outputs:
+        real ts_prof
+! local:
+        integer it_sandu1, it_sandu2,k
+        real timeit,time_sandu1,time_sandu2,frac
+! Check that initial day of the simulation consistent with SANDU period:
+       if (annee_ref.ne.2006 ) then
+        print*,'Pour SANDUREF, annee_ref doit etre 2006 '
+        print*,'Changer annee_ref dans run.def'
+        stop
+       endif
+!      if (annee_ref.eq.2006 .and. day1.lt.day_ini_sandu) then
+!       print*,'SANDUREF debute le 15 Juillet 2006 (jour julien=196)'
+!       print*,'Changer dayref dans run.def'
+!       stop
+!      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_sandu)*86400.
+!       else
+!        timeit=(day+61.-1.)*86400. ! 61 days between Nov01 and Dec31 1992
+!       endif
+      timeit=(day-day_ini_sandu)*86400
+
+! Determine the closest observation times:
+       it_sandu1=INT(timeit/dt_sandu)+1
+       it_sandu2=it_sandu1 + 1
+       time_sandu1=(it_sandu1-1)*dt_sandu
+       time_sandu2=(it_sandu2-1)*dt_sandu
+       print *,'timeit day day_ini_sandu',timeit,day,day_ini_sandu
+       print *,'it_sandu1,it_sandu2,time_sandu1,time_sandu2',
+     .          it_sandu1,it_sandu2,time_sandu1,time_sandu2
+
+       if (it_sandu1 .ge. nt_sandu) then
+        write(*,*) 'PB-stop: day, it_sandu1, it_sandu2, timeit: '
+     :        ,day,it_sandu1,it_sandu2,timeit/86400.
+        stop
+       endif
+
+! time interpolation:
+       frac=(time_sandu2-timeit)/(time_sandu2-time_sandu1)
+       frac=max(frac,0.0)
+
+       ts_prof = ts_sandu(it_sandu2)
+     :          -frac*(ts_sandu(it_sandu2)-ts_sandu(it_sandu1))
+
+         print*,
+     :'day,annee_ref,day_ini_sandu,timeit,it_sandu1,it_sandu2,SST:',
+     :day,annee_ref,day_ini_sandu,timeit/86400.,it_sandu1,
+     :it_sandu2,ts_prof
+
+        return
+        END
 !=====================================================================
 c-------------------------------------------------------------------------
@@ -2209 +2656 @@
           end
  
+!======================================================================
+        SUBROUTINE interp_astex_time(day,day1,annee_ref
+     i             ,year_ini_astex,day_ini_astex,nt_astex,dt_astex
+     i             ,nlev_astex,div_astex,ts_astex,ug_astex,vg_astex
+     i             ,ufa_astex,vfa_astex,div_prof,ts_prof,ug_prof,vg_prof
+     i             ,ufa_prof,vfa_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_astex: total nb of data in the forcing (e.g. 41 for Astex)
+! dt_astex: total time interval (in sec) between 2 forcing data (e.g. 1h for Astex)
+!---------------------------------------------------------------------------------------
+
+! inputs:
+        integer annee_ref
+        integer nt_astex,nlev_astex
+        integer year_ini_astex
+        real day, day1,day_ini_astex,dt_astex
+        real div_astex(nt_astex),ts_astex(nt_astex),ug_astex(nt_astex)
+        real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
+! outputs:
+        real div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
+! local:
+        integer it_astex1, it_astex2,k
+        real timeit,time_astex1,time_astex2,frac
+
+! Check that initial day of the simulation consistent with ASTEX period:
+       if (annee_ref.ne.1992 ) then
+        print*,'Pour Astex, annee_ref doit etre 1992 '
+        print*,'Changer annee_ref dans run.def'
+        stop
+       endif
+       if (annee_ref.eq.1992 .and. day1.lt.day_ini_astex) then
+        print*,'Astex debute le 13 Juin 1992 (jour julien=165)'
+        print*,'Changer dayref dans run.def'
+        stop
+       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_astex)*86400.
+!       else
+!        timeit=(day+2.-1.)*86400. ! 2 days between Jun13 and Jun15 1992
+!       endif
+      timeit=(day-day_ini_astex)*86400
+
+! Determine the closest observation times:
+       it_astex1=INT(timeit/dt_astex)+1
+       it_astex2=it_astex1 + 1
+       time_astex1=(it_astex1-1)*dt_astex
+       time_astex2=(it_astex2-1)*dt_astex
+       print *,'timeit day day_ini_astex',timeit,day,day_ini_astex
+       print *,'it_astex1,it_astex2,time_astex1,time_astex2',
+     .          it_astex1,it_astex2,time_astex1,time_astex2
+
+       if (it_astex1 .ge. nt_astex) then
+        write(*,*) 'PB-stop: day, it_astex1, it_astex2, timeit: '
+     :        ,day,it_astex1,it_astex2,timeit/86400.
+        stop
+       endif
+
+! time interpolation:
+       frac=(time_astex2-timeit)/(time_astex2-time_astex1)
+       frac=max(frac,0.0)
+
+       div_prof = div_astex(it_astex2)
+     :          -frac*(div_astex(it_astex2)-div_astex(it_astex1))
+       ts_prof = ts_astex(it_astex2)
+     :          -frac*(ts_astex(it_astex2)-ts_astex(it_astex1))
+       ug_prof = ug_astex(it_astex2)
+     :          -frac*(ug_astex(it_astex2)-ug_astex(it_astex1))
+       vg_prof = vg_astex(it_astex2)
+     :          -frac*(vg_astex(it_astex2)-vg_astex(it_astex1))
+       ufa_prof = ufa_astex(it_astex2)
+     :          -frac*(ufa_astex(it_astex2)-ufa_astex(it_astex1))
+       vfa_prof = vfa_astex(it_astex2)
+     :          -frac*(vfa_astex(it_astex2)-vfa_astex(it_astex1))
+
+         print*,
+     :'day,annee_ref,day_ini_astex,timeit,it_astex1,it_astex2,SST:',
+     :day,annee_ref,day_ini_astex,timeit/86400.,it_astex1,
+     :it_astex2,div_prof,ts_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
+
+        return
+        END
+
 !======================================================================
         SUBROUTINE interp_toga_time(day,day1,annee_ref
@@ -2479 +3017 @@
         return
         end
+!======================================================================
+      subroutine readprofile_sandu(nlev_max,kmax,height,pprof,tprof,
+     .       thlprof,qprof,uprof,vprof,wprof,omega,o3mmr)
+!======================================================================
+      implicit none
+
+        integer nlev_max,kmax,kmax2
+        logical :: llesread = .true.
+
+        real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
+     .  thlprof(nlev_max),
+     .  qprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
+     .  wprof(nlev_max),omega(nlev_max),o3mmr(nlev_max)
+
+        integer, parameter :: ilesfile=1
+        integer :: ierr,irad,imax,jtot,k
+        logical :: lmoist,lcoriol,ltimedep
+        real :: xsize,ysize
+        real :: ustin,wsvsurf,timerad
+        character(80) :: chmess
+
+        if(.not.(llesread)) return
+
+       open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
+        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
+        read (ilesfile,*) kmax
+        do k=1,kmax
+          read (ilesfile,*) height(k),pprof(k),  tprof(k),thlprof(k),
+     .                      qprof (k),uprof(k),  vprof(k),  wprof(k),
+     .                      omega (k),o3mmr(k)
+        enddo
+        close(ilesfile)
+
+        return
+        end
+
+!======================================================================
+      subroutine readprofile_astex(nlev_max,kmax,height,pprof,tprof,
+     .    thlprof,qvprof,qlprof,qtprof,uprof,vprof,wprof,tkeprof,o3mmr)
+!======================================================================
+      implicit none
+
+        integer nlev_max,kmax,kmax2
+        logical :: llesread = .true.
+
+        real height(nlev_max),pprof(nlev_max),tprof(nlev_max),
+     .  thlprof(nlev_max),qlprof(nlev_max),qtprof(nlev_max),
+     .  qvprof(nlev_max),uprof(nlev_max),vprof(nlev_max),
+     .  wprof(nlev_max),tkeprof(nlev_max),o3mmr(nlev_max)
+
+        integer, parameter :: ilesfile=1
+        integer :: ierr,irad,imax,jtot,k
+        logical :: lmoist,lcoriol,ltimedep
+        real :: xsize,ysize
+        real :: ustin,wsvsurf,timerad
+        character(80) :: chmess
+
+        if(.not.(llesread)) return
+
+       open (ilesfile,file='prof.inp.001',status='old',iostat=ierr)
+        if (ierr /= 0) stop 'ERROR:Prof.inp does not exist'
+        read (ilesfile,*) kmax
+        do k=1,kmax
+          read (ilesfile,*) height(k),pprof(k),  tprof(k),thlprof(k),
+     .                qvprof (k),qlprof (k),qtprof (k),
+     .                uprof(k),  vprof(k),  wprof(k),tkeprof(k),o3mmr(k)
+        enddo
+        close(ilesfile)
+
+        return
+        end
+
 
 
@@ -2539 +3149 @@
         return
         end
-!===============================================================
-      function ismin(n,sx,incx)
+!=====================================================================
+      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
-      integer n,i,incx,ismin,ix
-      real sx((n-1)*incx+1),sxmin
-
-      ix=1
-      ismin=1
-      sxmin=sx(1)
-      do i=1,n-1
-         ix=ix+incx
-         if(sx(ix).lt.sxmin) then
-             sxmin=sx(ix)
-             ismin=i+1
-         endif
-      enddo
-
-      return
-      end
-
-!===============================================================
-      function ismax(n,sx,incx)
+
+#include "netcdf.inc"
+
+      integer ntime,nlevel
+      integer l,k
+      character*80 :: fich_amma
+      real*8 time(ntime)
+      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
+     i         ,year_ini_amma,day_ini_amma,nt_amma,dt_amma,nlev_amma
+     i         ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma
+     o         ,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
+
+!=====================================================================
+      subroutine read_fire(fich_fire,nlevel,ntime
+     :     ,zz,thl,qt,u,v,tke
+     :     ,ug,vg,wls,dqtdx,dqtdy,dqtdt,thl_rad)
+
+!program reading forcings of the FIRE case study
+
 
       implicit none
-      integer n,i,incx,ismax,ix
-      real sx((n-1)*incx+1),sxmax
-
-      ix=1
-      ismax=1
-      sxmax=sx(1)
-      do i=1,n-1
-         ix=ix+incx
-         if(sx(ix).gt.sxmax) then
-             sxmax=sx(ix)
-             ismax=i+1
-         endif
-      enddo
-
-      return
-      end
-
+
+#include "netcdf.inc"
+
+      integer ntime,nlevel
+      integer l,k
+      character*80 :: fich_fire
+      real*8 time(ntime)
+      real*8 zz(nlevel) 
+
+      real*8 thl(nlevel)
+      real*8 qt(nlevel),u(nlevel)
+      real*8 v(nlevel),tke(nlevel)
+      real*8 ug(nlevel,ntime),vg(nlevel,ntime),wls(nlevel,ntime)
+      real*8 dqtdx(nlevel,ntime),dqtdy(nlevel,ntime)
+      real*8 dqtdt(nlevel,ntime),thl_rad(nlevel,ntime)  
+
+      integer nid, ierr
+      integer nbvar3d
+      parameter(nbvar3d=30)
+      integer var3didin(nbvar3d)
+
+      ierr = NF_OPEN(fich_fire,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,"thetal",var3didin(2))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'temp'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"qt",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,"tke",var3didin(6))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'tke'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"ugeo",var3didin(7))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'ug'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"vgeo",var3didin(8))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'vg'
+         endif
+      
+      ierr=NF_INQ_VARID(nid,"wls",var3didin(9))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'wls'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"dqtdx",var3didin(10))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'dqtdx'
+         endif
+
+      ierr=NF_INQ_VARID(nid,"dqtdy",var3didin(11))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'dqtdy'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"dqtdt",var3didin(12))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'dqtdt'
+      endif
+
+      ierr=NF_INQ_VARID(nid,"thl_rad",var3didin(13))
+         if(ierr/=NF_NOERR) then
+           write(*,*) NF_STRERROR(ierr)
+           stop 'thl_rad'
+      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),thl)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(2),thl)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture thl ok',thl
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(3),qt)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(3),qt)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture qt ok',qt
+ 
+#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),tke)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(6),tke)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture tke ok',tke
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(7),ug)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(7),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(8),vg)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(8),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(9),wls)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(9),wls)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture wls ok',wls
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(10),dqtdx)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(10),dqtdx)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture dqtdx ok',dqtdx
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(11),dqtdy)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(11),dqtdy)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture dqtdy ok',dqtdy
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(12),dqtdt)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(12),dqtdt)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture dqtdt ok',dqtdt
+
+#ifdef NC_DOUBLE
+         ierr = NF_GET_VAR_DOUBLE(nid,var3didin(13),thl_rad)
+#else
+         ierr = NF_GET_VAR_REAL(nid,var3didin(13),thl_rad)
+#endif
+         if(ierr/=NF_NOERR) then
+            write(*,*) NF_STRERROR(ierr)
+            stop "getvarup"
+         endif
+!          write(*,*)'lecture thl_rad ok',thl_rad
+
+         return 
+         end subroutine read_fire