Changeset 1795 for LMDZ5/branches/testing/libf/phy1d

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)
• libf/phy1d/1DUTILS.h_with_writelim (modified) (4 diffs)
• libf/phy1d/1DUTILS.h_with_writelim_old (modified) (4 diffs)
• libf/phy1d/1D_decl_cases.h (modified) (2 diffs)
• libf/phy1d/1D_interp_cases.h (modified) (2 diffs)
• libf/phy1d/1D_read_forc_cases.h (modified) (8 diffs)
• libf/phy1d/compar1d.h (modified) (1 diff)
• libf/phy1d/lmdz1d.F (modified) (21 diffs)
• libf/phy1d/ocean_forced_mod.F90 (modified) (2 diffs)
• libf/phy1d/surf_land_bucket_mod.F90 (modified) (1 diff)
• libf/phy1d/tsoilnudge.h (copied) (copied from LMDZ5/trunk/libf/phy1d/tsoilnudge.h)

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

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

@@ -125 +125 @@
        ok_flux_surf = .FALSE.
        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
@@ -703 +710 @@
       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
@@ -1079 +1069 @@
  
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-      SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
+      SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
  
 !    Auteur :  P. Le Van .
@@ -2662 +2652 @@
         return
         end
-!===============================================================
-      function ismin(n,sx,incx)
-
-      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)
-
-      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
-
+

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

@@ -125 +125 @@
        ok_flux_surf = .FALSE.
        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
@@ -703 +710 @@
       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
@@ -1079 +1069 @@
  
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-      SUBROUTINE disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
+      SUBROUTINE disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
  
 !    Auteur :  P. Le Van .
@@ -2604 +2594 @@
         end
 
-!===============================================================
-      function ismin(n,sx,incx)
-
-      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)
-
-      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
-
+

LMDZ5/branches/testing/libf/phy1d/1D_decl_cases.h

@@ -81 +81 @@
         real ht_proftwp(nlev_twpi),vt_proftwp(nlev_twpi)
         real hq_proftwp(nlev_twpi),vq_proftwp(nlev_twpi)
+
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!Declarations specifiques au cas AMMA
+        character*80 :: fich_amma
+! Option du cas AMMA ou on impose la discretisation verticale (Ap,Bp)
+        logical  :: fixe_disvert=.true.
+        integer nlev_amma, nt_amma
+!       parameter (nlev_amma=29, nt_amma=48)  ! Fleur, juillet 2012
+        parameter (nlev_amma=36, nt_amma=48)  ! Romain, octobre 2012
+!       parameter (nlev_amma=26, nt_amma=48)  ! Test MPL feverier 2013
+        integer year_ini_amma, day_ini_amma, mth_ini_amma
+        real heure_ini_amma
+        real day_ju_ini_amma   ! Julian day of amma first day
+        parameter (year_ini_amma=2006) 
+        parameter (mth_ini_amma=7)
+        parameter (day_ini_amma=10)  ! 10 = 10Juil2006
+        parameter (heure_ini_amma=0.) !0h en secondes
+        real dt_amma
+        parameter (dt_amma=1800.)
+
+!profils initiaux:
+        real plev_amma(nlev_amma)
+        real tv_amma(nlev_amma),rho_amma(nlev_amma)
+        real thv_amma(nlev_amma)
+        
+        real z_amma(nlev_amma)
+        real th_amma(nlev_amma),q_amma(nlev_amma)
+        real u_amma(nlev_amma)
+        real v_amma(nlev_amma)
+
+        real thvsurf_amma,tvsurf_amma,rhosurf_amma,thsurf
+
+        real th_ammai(nlev_amma),q_ammai(nlev_amma)
+        real u_ammai(nlev_amma)
+        real v_ammai(nlev_amma)
+        real vitw_ammai(nlev_amma) 
+        real ht_ammai(nlev_amma)
+        real hq_ammai(nlev_amma)
+        real vt_ammai(nlev_amma)
+        real vq_ammai(nlev_amma)
+        
+!forcings
+        real ht_amma(nlev_amma,nt_amma)
+        real hq_amma(nlev_amma,nt_amma)
+        real vitw_amma(nlev_amma,nt_amma)
+        real lat_amma(nt_amma),sens_amma(nt_amma)
+
+!champs interpoles
+        real plev_profamma(nlev_amma),vitw_profamma(nlev_amma)
+        real ht_profamma(nlev_amma)
+        real hq_profamma(nlev_amma)
+        real lat_profamma,sens_profamma
+        real vt_profamma(nlev_amma)
+        real vq_profamma(nlev_amma)
+        real th_profamma(nlev_amma)
+        real q_profamma(nlev_amma)
+        real u_profamma(nlev_amma)
+        real v_profamma(nlev_amma)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!Declarations specifiques au cas FIRE
+        character*80 :: fich_fire
+        integer nlev_fire, nt_fire
+        parameter (nlev_fire=120, nt_fire=1)  
+        integer year_ini_fire, day_ini_fire, mth_ini_fire
+        real heure_ini_fire
+        real day_ju_ini_fire   ! Julian day of fire first day
+        parameter (year_ini_fire=1987) 
+        parameter (mth_ini_fire=7)
+        parameter (day_ini_fire=14)  ! 14 = 14Juil1987
+        parameter (heure_ini_fire=0.) !0h en secondes
+
+!profils initiaux:
+        real z_fire(nlev_fire)
+        real thl_fire(nlev_fire),qt_fire(nlev_fire)
+        real u_fire(nlev_fire), v_fire(nlev_fire)
+        real tke_fire(nlev_fire) 
+        
+!forcings
+        real ugeo_fire(nlev_fire),vgeo_fire(nlev_fire)
+        real wls_fire(nlev_fire),dqtdx_fire(nlev_fire)
+        real dqtdy_fire(nlev_fire)
+        real dqtdt_fire(nlev_fire),thl_rad_fire(nlev_fire)
+         
+        real ugeo_mod(llm),vgeo_mod(llm),wls_mod(llm)
+        real dqtdx_mod(llm),dqtdy_mod(llm),dqtdt_mod(llm)
+        real thl_rad_mod(llm)
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! Declarations specifiques au cas GCSSold
@@ -127 +215 @@
         real sens_prof,flat_prof,fact
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+! declarations specifiques au cas Sandu
+        character*80 :: fich_sandu
+!        integer nlev_prof
+!        parameter (nlev_prof = 41)
+        integer nlev_sandu, nt_sandu
+        parameter (nlev_sandu=87, nt_sandu=13)
+        integer year_ini_sandu, day_ini_sandu, mth_ini_sandu
+        real day_ju_ini_sandu                                ! Julian day of sandu case first day
+        parameter (year_ini_sandu=2006)
+        parameter (mth_ini_sandu=7)
+        parameter (day_ini_sandu=15)  ! 196 = 15 juillet 2006
+        real dt_sandu, tau_sandu
+        logical  :: trouve_700=.true.
+        parameter (dt_sandu=6.*3600.)   ! forcages donnes ttes les 6 heures par ifa_sandu.txt
+!       parameter (tau_sandu=3600.)  ! temps de relaxation u,v,thetal,qt vers profil init et au dessus 700hPa
+!!
+        integer it_sandu1, it_sandu2
+        real time_sandu1,time_sandu2
+
+        real ts_sandu(nt_sandu)
+! profs comme "profil sandu"
+        real plev_profs(nlev_sandu)
+        real t_profs(nlev_sandu),thl_profs(nlev_sandu)
+        real q_profs(nlev_sandu)
+        real u_profs(nlev_sandu),v_profs(nlev_sandu),w_profs(nlev_sandu)
+        real omega_profs(nlev_sandu),o3mmr_profs(nlev_sandu)
+
+        real thl_mod(llm),omega_mod(llm),o3mmr_mod(llm),tke_mod(llm)
+! pour relaxer u,v,thl et qt vers les profils initiaux au dessus de 700hPa
+        real relax_u(llm),relax_v(llm),relax_thl(llm),relax_q(llm,2)
+!vertical advection computation
+        real d_t_z(llm), d_q_z(llm)
+        real d_t_dyn_z(llm), d_q_dyn_z(llm)
+        real zz(llm)
+        real zfact
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Declarations specifiques au cas Astex
+        character*80 :: fich_astex
+        integer nlev_astex, nt_astex
+        parameter (nlev_astex=34, nt_astex=49)
+        integer year_ini_astex, day_ini_astex, mth_ini_astex
+        real day_ju_ini_astex                                ! Julian day of astex case first day
+        parameter (year_ini_astex=1992)
+        parameter (mth_ini_astex=6)
+        parameter (day_ini_astex=13)  ! 165 = 13 juin 1992
+        real dt_astex, tau_astex
+        parameter (dt_astex=3600.)    ! forcages donnes ttes les heures par ifa_astex.txt
+        integer it_astex1, it_astex2
+        real time_astex1,time_astex2
+        real ts_astex(nt_astex),div_astex(nt_astex),ug_astex(nt_astex)
+        real vg_astex(nt_astex),ufa_astex(nt_astex),vfa_astex(nt_astex)
+        real div_prof,ug_prof,vg_prof,ufa_prof,vfa_prof
+! profa comme "profil astex"
+        real plev_profa(nlev_astex)
+        real t_profa(nlev_astex),thl_profa(nlev_astex)
+        real qv_profa(nlev_astex),ql_profa(nlev_astex)
+        real qt_profa(nlev_astex),o3mmr_profa(nlev_astex)
+        real u_profa(nlev_astex),v_profa(nlev_astex),w_profa(nlev_astex)
+        real tke_profa(nlev_astex)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+

LMDZ5/branches/testing/libf/phy1d/1D_interp_cases.h

@@ -181 +181 @@
 
       endif ! forcing_twpice
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
+! Interpolation forcing AMMA
+!---------------------------------------------------------------------
+
+       if (forcing_amma) then
+
+        print*,
+     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_amma=',
+     :    daytime,day1,(daytime-day1)*86400.,
+     :    (daytime-day1)*86400/dt_amma
+
+! time interpolation using TOGA interpolation routine
+        CALL interp_amma_time(daytime,day1,annee_ref
+     i       ,year_ini_amma,day_ju_ini_amma,nt_amma,dt_amma,nlev_amma
+     i       ,vitw_amma,ht_amma,hq_amma,lat_amma,sens_amma
+     o       ,vitw_profamma,ht_profamma,hq_profamma,lat_profamma
+     :       ,sens_profamma)
+
+      print*,'apres interpolation temporelle AMMA'
+
+      do k=1,nlev_amma
+         th_profamma(k)=0.
+         q_profamma(k)=0.
+         u_profamma(k)=0.
+         v_profamma(k)=0.
+         vt_profamma(k)=0.
+         vq_profamma(k)=0.
+       enddo
+! vertical interpolation using TOGA interpolation routine:
+!      write(*,*)'avant interp vert', t_proftwp
+      CALL interp_toga_vertical(play,nlev_amma,plev_amma
+     :         ,th_profamma,q_profamma,u_profamma,v_profamma
+     :         ,vitw_profamma
+     :         ,ht_profamma,vt_profamma,hq_profamma,vq_profamma
+     :         ,t_mod,q_mod,u_mod,v_mod,w_mod
+     :         ,ht_mod,vt_mod,hq_mod,vq_mod,mxcalc)
+       write(*,*) 'Profil initial forcing AMMA interpole'
+
+
+!calcul de l'advection verticale a partir du omega
+cCalcul des gradients verticaux
+cinitialisation
+      do l=1,llm
+      d_t_z(l)=0.
+      d_q_z(l)=0.
+      enddo
+
+      DO l=2,llm-1
+       d_t_z(l)=(temp(l+1)-temp(l-1))
+     &          /(play(l+1)-play(l-1))
+       d_q_z(l)=(q(l+1,1)-q(l-1,1))
+     &          /(play(l+1)-play(l-1))
+      ENDDO
+      d_t_z(1)=d_t_z(2)
+      d_q_z(1)=d_q_z(2)
+      d_t_z(llm)=d_t_z(llm-1)
+      d_q_z(llm)=d_q_z(llm-1)
+
+
+      do l = 1, llm
+       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
+       omega(l) = w_mod(l)*(-rg*rho(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)
+!calcul de l'advection totale
+!        d_th_adv(l) = alpha*omega(l)/rcpd+ht_mod(l)-omega(l)*d_t_z(l)
+!attention: on impose dth
+        d_th_adv(l) = alpha*omega(l)/rcpd+
+     &         ht_mod(l)*(play(l)/pzero)**rkappa-omega(l)*d_t_z(l)
+!        d_th_adv(l) = 0.
+!        print*,'temp vert adv',l,ht_mod(l),vt_mod(l),-d_t_dyn_z(l)
+        d_q_adv(l,1) = hq_mod(l)-omega(l)*d_q_z(l)
+!        d_q_adv(l,1) = 0.
+!        print*,'q vert adv',l,hq_mod(l),vq_mod(l),-d_q_dyn_z(l)
+   
+       dt_cooling(l) = 0.0
+      enddo
+
+
+!     ok_flux_surf=.false.
+      fsens=-1.*sens_profamma
+      flat=-1.*lat_profamma
+
+      endif ! forcing_amma
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !---------------------------------------------------------------------
@@ -254 +341 @@
       endif ! forcing_armcu
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+!---------------------------------------------------------------------
+! Interpolation forcing in time and onto model levels
+!---------------------------------------------------------------------
+      if (forcing_sandu) then
+
+        print*,
+     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_sandu=',
+     :    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_sandu
+
+! time interpolation:
+! ATTENTION, cet appel ne convient pas pour TOGA !!
+! revoir 1DUTILS.h et les arguments
+      CALL interp_sandu_time(daytime,day1,annee_ref
+     i             ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu
+     i             ,nlev_sandu
+     i             ,ts_sandu,ts_prof)
+
+        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
+
+! vertical interpolation:
+      CALL interp_sandu_vertical(play,nlev_sandu,plev_profs
+     :         ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs
+     :         ,omega_profs,o3mmr_profs
+     :         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod
+     :         ,omega_mod,o3mmr_mod,mxcalc)
+!calcul de l'advection verticale
+cCalcul des gradients verticaux
+cinitialisation
+      d_t_z(:)=0.
+      d_q_z(:)=0.
+      d_t_dyn_z(:)=0.
+      d_q_dyn_z(:)=0.
+! schema centre
+!     DO l=2,llm-1
+!      d_t_z(l)=(temp(l+1)-temp(l-1))
+!    &          /(play(l+1)-play(l-1))
+!      d_q_z(l)=(q(l+1,1)-q(l-1,1))
+!    &          /(play(l+1)-play(l-1))
+! schema amont
+      DO l=2,llm-1
+       d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
+       d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
+!     print *,'l temp2 temp0 play2 play0 omega_mod',
+!    & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l)
+      ENDDO
+      d_t_z(1)=d_t_z(2)
+      d_q_z(1)=d_q_z(2)
+      d_t_z(llm)=d_t_z(llm-1)
+      d_q_z(llm)=d_q_z(llm-1)
+
+!  calcul de l advection verticale
+! Confusion w (m/s) et omega (Pa/s) !!
+      d_t_dyn_z(:)=omega_mod(:)*d_t_z(:)
+      d_q_dyn_z(:)=omega_mod(:)*d_q_z(:)
+!     do l=1,llm
+!      print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
+!    :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l)
+!     enddo
+
+
+! large-scale forcing : pour le cas Sandu ces forcages sont la SST
+! et une divergence constante -> profil de omega
+      tsurf = ts_prof
+      write(*,*) 'SST suivante: ',tsurf
+      do l = 1, llm
+       omega(l) = omega_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)
+!
+!      d_th_adv(l) = 0.0
+!      d_q_adv(l,1) = 0.0
+!CR:test advection=0
+!calcul de l'advection verticale
+        d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
+!        print*,'temp adv',l,-d_t_dyn_z(l)
+        d_q_adv(l,1) = -d_q_dyn_z(l)
+!        print*,'q adv',l,-d_q_dyn_z(l)
+       dt_cooling(l) = 0.0
+      enddo
+      endif ! forcing_sandu
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
+! Interpolation forcing in time and onto model levels
+!---------------------------------------------------------------------
+      if (forcing_astex) then
+
+        print*,
+     : '#### ITAP,day,day1,(day-day1)*86400,(day-day1)*86400/dt_astex=',
+     :    day,day1,(day-day1)*86400.,(day-day1)*86400/dt_astex
+
+! time interpolation:
+! ATTENTION, cet appel ne convient pas pour TOGA !!
+! revoir 1DUTILS.h et les arguments
+      CALL interp_astex_time(daytime,day1,annee_ref
+     i             ,year_ini_astex,day_ju_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)
+
+        if (type_ts_forcing.eq.1) ts_cur = ts_prof ! SST used in read_tsurf1d
+
+! vertical interpolation:
+      CALL interp_astex_vertical(play,nlev_astex,plev_profa
+     :         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa
+     :         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa
+     :         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod
+     :         ,tke_mod,o3mmr_mod,mxcalc)
+!calcul de l'advection verticale
+!Calcul des gradients verticaux
+!initialisation
+      d_t_z(:)=0.
+      d_q_z(:)=0.
+      d_t_dyn_z(:)=0.
+      d_q_dyn_z(:)=0.
+! schema centre
+!     DO l=2,llm-1
+!      d_t_z(l)=(temp(l+1)-temp(l-1))
+!    &          /(play(l+1)-play(l-1))
+!      d_q_z(l)=(q(l+1,1)-q(l-1,1))
+!    &          /(play(l+1)-play(l-1))
+! schema amont
+      DO l=2,llm-1
+       d_t_z(l)=(temp(l+1)-temp(l))/(play(l+1)-play(l))
+       d_q_z(l)=(q(l+1,1)-q(l,1))/(play(l+1)-play(l))
+!     print *,'l temp2 temp0 play2 play0 omega_mod',
+!    & temp(l+1),temp(l-1),play(l+1),play(l-1),omega_mod(l)
+      ENDDO
+      d_t_z(1)=d_t_z(2)
+      d_q_z(1)=d_q_z(2)
+      d_t_z(llm)=d_t_z(llm-1)
+      d_q_z(llm)=d_q_z(llm-1)
+
+!  calcul de l advection verticale
+! Confusion w (m/s) et omega (Pa/s) !!
+      d_t_dyn_z(:)=w_mod(:)*d_t_z(:)
+      d_q_dyn_z(:)=w_mod(:)*d_q_z(:)
+!     do l=1,llm
+!      print *,'d_t_dyn omega_mod d_t_z d_q_dyn d_q_z',
+!    :l,d_t_dyn_z(l),omega_mod(l),d_t_z(l),d_q_dyn_z(l),d_q_z(l)
+!     enddo
+
+
+! large-scale forcing : pour le cas Astex ces forcages sont la SST
+! la divergence,ug,vg,ufa,vfa
+      tsurf = ts_prof
+      write(*,*) 'SST suivante: ',tsurf
+      do l = 1, llm
+       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)
+!
+!      d_th_adv(l) = 0.0
+!      d_q_adv(l,1) = 0.0
+!CR:test advection=0
+!calcul de l'advection verticale
+        d_th_adv(l) = alpha*omega(l)/rcpd-d_t_dyn_z(l)
+!        print*,'temp adv',l,-d_t_dyn_z(l)
+        d_q_adv(l,1) = -d_q_dyn_z(l)
+!        print*,'q adv',l,-d_q_dyn_z(l)
+       dt_cooling(l) = 0.0
+      enddo
+      endif ! forcing_astex
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+

LMDZ5/branches/testing/libf/phy1d/1D_read_forc_cases.h

@@ -4 +4 @@
 !----------------------------------------------------------------------
 
-      if (forcing_les .or. forcing_radconv .or. forcing_GCSSold ) then
-
+      if (forcing_les .or. forcing_radconv
+     :    .or. forcing_GCSSold .or. forcing_fire) then
+
+      if (forcing_fire) then
+!----------------------------------------------------------------------
+!read fire forcings from fire.nc
+!----------------------------------------------------------------------
+      fich_fire='fire.nc'
+      call read_fire(fich_fire,nlev_fire,nt_fire
+     :     ,height,tttprof,qtprof,uprof,vprof,e12prof
+     :     ,ugprof,vgprof,wfls,dqtdxls
+     :     ,dqtdyls,dqtdtls,thlpcar)
+      write(*,*) 'Forcing FIRE lu'
+      kmax=120            ! nombre de niveaux dans les profils et forcages
+      else 
 !----------------------------------------------------------------------
 ! Read profiles from files: prof.inp.001 and lscale.inp.001
@@ -16 +29 @@
      .           wfls,dqtdxls,dqtdyls,dqtdtls,
      .           thlpcar) 
+      endif
 
 ! compute altitudes of play levels.
@@ -35 +49 @@
         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
+       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
           temp(l) = ttt*(play(l)/pzero)**rkappa
           teta(l) = ttt
-        else
+       else
           temp(l) = ttt
           teta(l) = ttt*(pzero/play(l))**rkappa
-        endif
+       endif
           print *,' temp,teta ',l,temp(l),teta(l)
         q(l,1)  = qtprof(kmax)-frac*( qtprof(kmax)- qtprof(kmax-1))
@@ -58 +72 @@
           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
+       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
           temp(l) = ttt*(play(l)/pzero)**rkappa
           teta(l) = ttt
-        else
+       else
           temp(l) = ttt
           teta(l) = ttt*(pzero/play(l))**rkappa
-        endif
+       endif
           print *,' temp,teta ',l,temp(l),teta(l)
             q(l,1)  = qtprof(k)-frac*( qtprof(k)- qtprof(k-1))
@@ -77 +91 @@
           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
+       if ((forcing_GCSSold .AND. tp_ini_GCSSold) .OR. forcing_fire)then ! pot. temp. in initial profile
           temp(l) = ttt*(play(l)/pzero)**rkappa
           teta(l) = ttt
-        else
+       else
           temp(l) = ttt
           teta(l) = ttt*(pzero/play(l))**rkappa
-        endif
+       endif
             q(l,1)  = qtprof(1)
             u(l)    =  uprof(1)
@@ -112 +126 @@
       enddo
 
-      endif ! forcing_les .or. forcing_GCSSold
+      endif ! forcing_les .or. forcing_GCSSold .or. forcing_fire
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !---------------------------------------------------------------------
@@ -263 +277 @@
        
       endif !forcing_twpice
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
+! Forcing from AMMA experiment (Couvreux et al. 2010) :
+!---------------------------------------------------------------------
+
+      if (forcing_amma) then
+!read AMMA forcings
+      fich_amma='amma.nc'
+      call read_amma(fich_amma,nlev_amma,nt_amma
+     :     ,z_amma,plev_amma,th_amma,q_amma,u_amma,v_amma,vitw_amma
+     :     ,ht_amma,hq_amma,sens_amma,lat_amma)
+
+      write(*,*) 'Forcing AMMA lu'
+
+!champs initiaux:
+      do k=1,nlev_amma
+         th_ammai(k)=th_amma(k)
+         q_ammai(k)=q_amma(k)
+         u_ammai(k)=u_amma(k)
+         v_ammai(k)=v_amma(k)
+         vitw_ammai(k)=vitw_amma(k,12)
+         ht_ammai(k)=ht_amma(k,12)
+         hq_ammai(k)=hq_amma(k,12)
+         vt_ammai(k)=0.
+         vq_ammai(k)=0.
+      enddo   
+      omega(:)=0.      
+      omega2(:)=0.
+      rho(:)=0.
+! vertical interpolation using TOGA interpolation routine:
+!      write(*,*)'avant interp vert', t_proftwp
+      CALL interp_toga_vertical(play,nlev_amma,plev_amma
+     :         ,th_ammai,q_ammai,u_ammai,v_ammai,vitw_ammai
+     :         ,ht_ammai,vt_ammai,hq_ammai,vq_ammai
+     :         ,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 mxcalc: ',tsurf,mxcalc
+      do l = 1, llm
+! Ligne du dessous à decommenter si on lit theta au lieu de temp
+!      temp(l) = t_mod(l)*(play(l)/pzero)**rkappa 
+       temp(l) = t_mod(l) 
+       q(l,1) = q_mod(l)
+       q(l,2) = 0.0
+!      print *,'read_forc: l,temp,q=',l,temp(l),q(l,1)
+       u(l) = u_mod(l)
+       v(l) = v_mod(l)
+       rho(l)  = play(l)/(rd*temp(l)*(1.+(rv/rd-1.)*q(l,1)))
+       omega(l) = w_mod(l)*(-rg*rho(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)
+!forcage en th
+!       d_th_adv(l) = ht_mod(l)
+       d_q_adv(l,1) = hq_mod(l)
+       d_q_adv(l,2) = 0.0
+       dt_cooling(l)=0.
+      enddo     
+       write(*,*) 'Profil initial forcing AMMA interpole temp39',
+     &           temp(39)
+      
+
+!     ok_flux_surf=.false.
+      fsens=-1.*sens_amma(12)
+      flat=-1.*lat_amma(12)
+       
+      endif !forcing_amma
+
+
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !---------------------------------------------------------------------
@@ -366 +456 @@
       endif ! forcing_armcu
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
+!---------------------------------------------------------------------
+! Forcing from transition case of Irina Sandu                  
+!---------------------------------------------------------------------
+
+      if (forcing_sandu) then
+       write(*,*) 'Avant lecture Forcing SANDU'
+
+! read sanduref forcing :
+      fich_sandu = './ifa_sanduref.txt'
+      CALL read_sandu(fich_sandu,nlev_sandu,nt_sandu,ts_sandu)
+
+       write(*,*) 'Forcing SANDU lu'
+
+!----------------------------------------------------------------------
+! Read profiles from file: prof.inp.001 
+!----------------------------------------------------------------------
+
+      call readprofile_sandu(nlev_max,kmax,height,plev_profs,t_profs,
+     .           thl_profs,q_profs,u_profs,v_profs,
+     .           w_profs,omega_profs,o3mmr_profs)
+
+! time interpolation for initial conditions:
+      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
+! ATTENTION, cet appel ne convient pas pour le cas SANDU !!
+! revoir 1DUTILS.h et les arguments
+
+      print *,'Avant interp_sandu_time'
+      print *,'daytime=',daytime
+      print *,'day1=',day1
+      print *,'annee_ref=',annee_ref
+      print *,'year_ini_sandu=',year_ini_sandu
+      print *,'day_ju_ini_sandu=',day_ju_ini_sandu
+      print *,'nt_sandu=',nt_sandu
+      print *,'dt_sandu=',dt_sandu
+      print *,'nlev_sandu=',nlev_sandu
+      CALL interp_sandu_time(daytime,day1,annee_ref
+     i             ,year_ini_sandu,day_ju_ini_sandu,nt_sandu,dt_sandu
+     i             ,nlev_sandu
+     i             ,ts_sandu,ts_prof)
+
+! vertical interpolation:
+      print *,'Avant interp_vertical: nlev_sandu=',nlev_sandu
+      CALL interp_sandu_vertical(play,nlev_sandu,plev_profs
+     :         ,t_profs,thl_profs,q_profs,u_profs,v_profs,w_profs
+     :         ,omega_profs,o3mmr_profs
+     :         ,t_mod,thl_mod,q_mod,u_mod,v_mod,w_mod
+     :         ,omega_mod,o3mmr_mod,mxcalc)
+       write(*,*) 'Profil initial forcing SANDU interpole'
+
+! initial and boundary conditions :
+      tsurf = ts_prof
+      write(*,*) 'SST initiale: ',tsurf
+      do l = 1, llm
+       temp(l) = t_mod(l)
+       tetal(l)=thl_mod(l)
+       q(l,1) = q_mod(l)
+       q(l,2) = 0.0
+       u(l) = u_mod(l)
+       v(l) = v_mod(l)
+       w(l) = w_mod(l)
+       omega(l) = omega_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+vt_mod(l)
+!      d_q_adv(l,1) = vq_mod(l)
+       d_th_adv(l) = alpha*omega(l)/rcpd
+       d_q_adv(l,1) = 0.0
+       d_q_adv(l,2) = 0.0
+      enddo
+
+      endif ! forcing_sandu
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!---------------------------------------------------------------------
+! Forcing from Astex case
+!---------------------------------------------------------------------
+
+      if (forcing_astex) then
+       write(*,*) 'Avant lecture Forcing Astex'
+
+! read astex forcing :
+      fich_astex = './ifa_astex.txt'
+      CALL read_astex(fich_astex,nlev_astex,nt_astex,div_astex,ts_astex,
+     :  ug_astex,vg_astex,ufa_astex,vfa_astex)
+
+       write(*,*) 'Forcing Astex lu'
+
+!----------------------------------------------------------------------
+! Read profiles from file: prof.inp.001
+!----------------------------------------------------------------------
+
+      call readprofile_astex(nlev_max,kmax,height,plev_profa,t_profa,
+     .           thl_profa,qv_profa,ql_profa,qt_profa,u_profa,v_profa,
+     .           w_profa,tke_profa,o3mmr_profa)
+
+! time interpolation for initial conditions:
+      write(*,*) 'AVT 1ere INTERPOLATION: day,day1 = ',day,day1
+! ATTENTION, cet appel ne convient pas pour le cas SANDU !!
+! revoir 1DUTILS.h et les arguments
+
+      print *,'Avant interp_astex_time'
+      print *,'daytime=',daytime
+      print *,'day1=',day1
+      print *,'annee_ref=',annee_ref
+      print *,'year_ini_astex=',year_ini_astex
+      print *,'day_ju_ini_astex=',day_ju_ini_astex
+      print *,'nt_astex=',nt_astex
+      print *,'dt_astex=',dt_astex
+      print *,'nlev_astex=',nlev_astex
+      CALL interp_astex_time(daytime,day1,annee_ref
+     i             ,year_ini_astex,day_ju_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)
+
+! vertical interpolation:
+      print *,'Avant interp_vertical: nlev_astex=',nlev_astex
+      CALL interp_astex_vertical(play,nlev_astex,plev_profa
+     :         ,t_profa,thl_profa,qv_profa,ql_profa,qt_profa
+     :         ,u_profa,v_profa,w_profa,tke_profa,o3mmr_profa
+     :         ,t_mod,thl_mod,qv_mod,ql_mod,qt_mod,u_mod,v_mod,w_mod
+     :         ,tke_mod,o3mmr_mod,mxcalc)
+       write(*,*) 'Profil initial forcing Astex interpole'
+
+! initial and boundary conditions :
+      tsurf = ts_prof
+      write(*,*) 'SST initiale: ',tsurf
+      do l = 1, llm
+       temp(l) = t_mod(l)
+       tetal(l)=thl_mod(l)
+       q(l,1) = qv_mod(l)
+       q(l,2) = ql_mod(l)
+       u(l) = u_mod(l)
+       v(l) = v_mod(l)
+       w(l) = w_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+vt_mod(l)
+!      d_q_adv(l,1) = vq_mod(l)
+       d_th_adv(l) = alpha*omega(l)/rcpd
+       d_q_adv(l,1) = 0.0
+       d_q_adv(l,2) = 0.0
+      enddo
+
+      endif ! forcing_astex
+

LMDZ5/branches/testing/libf/phy1d/compar1d.h

@@ -25 +25 @@
 
       logical :: restart
+      logical :: ok_old_disvert
 
       common/com_par1d/forcing_type,nat_surf,tsurf,rugos,               &
      & qsol,qsurf,psurf,zsurf,albedo,time,time_ini,xlat,xlon,airefi,    &
      & wtsurf,wqsurf,restart_runoff,xagesno,qsolinp,zpicinp,            &
-     & restart
+     & restart,ok_old_disvert
 
 !$OMP THREADPRIVATE(/com_par1d/)

LMDZ5/branches/testing/libf/phy1d/lmdz1d.F

@@ -6 +6 @@
       use dimphy
       use surface_data, only : type_ocean,ok_veget
-      use pbl_surface_mod, only : pbl_surface_init, pbl_surface_final
+      use pbl_surface_mod, only : ftsoil, pbl_surface_init,
+     $                            pbl_surface_final
       use fonte_neige_mod, only : fonte_neige_init, fonte_neige_final 
 
       use infotrac ! new
       use control_mod
+      USE indice_sol_mod
 
       implicit none
@@ -20 +22 @@
 #include "clesphys.h"
 #include "dimsoil.h"
-#include "indicesol.h"
+!#include "indicesol.h"
 
 #include "comvert.h"
 #include "compar1d.h"
 #include "flux_arp.h"
+#include "tsoilnudge.h"
 #include "fcg_gcssold.h"
 !!!#include "fbforcing.h"
@@ -86 +89 @@
 
         integer :: kmax = llm
-        integer nlev_max
-        parameter (nlev_max = 100)
+        integer nlev_max,llm700
+        parameter (nlev_max = 1000)
         real timestep, frac, timeit
         real height(nlev_max),tttprof(nlev_max),qtprof(nlev_max),
@@ -98 +101 @@
 c        integer :: forcing_type
         logical :: forcing_les     = .false.
-        logical :: forcing_armcu  = .false.
+        logical :: forcing_armcu   = .false.
         logical :: forcing_rico    = .false.
         logical :: forcing_radconv = .false.
         logical :: forcing_toga    = .false.
         logical :: forcing_twpice  = .false.
+        logical :: forcing_amma    = .false.
         logical :: forcing_GCM2SCM = .false.
         logical :: forcing_GCSSold = .false.
+        logical :: forcing_sandu   = .false.
+        logical :: forcing_astex   = .false.
+        logical :: forcing_fire    = .false.
         integer :: type_ts_forcing ! 0 = SST constant; 1 = SST read from a file
 !                                                            (cf read_tsurf1d.F)
 
 !vertical advection computation
-        real d_t_z(llm), d_q_z(llm)
-        real d_t_dyn_z(llm), d_q_dyn_z(llm)
-        real zz(llm)
-        real zfact
+!       real d_t_z(llm), d_q_z(llm)
+!       real d_t_dyn_z(llm), d_q_dyn_z(llm)
+!       real zz(llm)
+!       real zfact
 
 !flag forcings
@@ -129 +136 @@
       real :: pzero=1.e5
       real :: play (llm),zlay (llm),sig_s(llm),plev(llm+1)
-      real :: playd(llm),zlayd(llm)
+      real :: playd(llm),zlayd(llm),ap_amma(llm+1),bp_amma(llm+1),poub
 
 !---------------------------------------------------------------------
@@ -137 +144 @@
       integer :: iq
       real :: phi(llm)
-      real :: teta(llm),temp(llm),u(llm),v(llm)
+      real :: teta(llm),tetal(llm),temp(llm),u(llm),v(llm),w(llm)
+      real :: rlat_rad(1),rlon_rad(1)
       real :: omega(llm+1),omega2(llm),rho(llm+1)
       real :: ug(llm),vg(llm),fcoriolis
+      real :: sfdt, cfdt
       real :: du_phys(llm),dv_phys(llm),dt_phys(llm)
       real :: du_dyn(llm),dv_dyn(llm),dt_dyn(llm)
@@ -194 +203 @@
 !  Fichiers et d'autres variables
 !---------------------------------------------------------------------
-      real ttt
+      real ttt,bow,q1
       integer :: ierr,k,l,i,it=1,mxcalc
       integer jjmp1
@@ -250 +259 @@
 !             initial profiles from RICO files 
 !             LS convergence imposed from RICO files 
+!forcing_type = 6 ==> forcing_amma = .true.
+!             initial profiles from AMMA nc file 
+!             LS convergence, omega and surface fluxes imposed from AMMA file  
 !forcing_type = 40 ==> forcing_GCSSold = .true.
 !             initial profile from GCSS file
 !             LS convergence imposed from GCSS file
+!forcing_type = 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
+!forcing_type = 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
 !forcing_type = 61 ==> forcing_armcu = .true.
-!             initial profile from arm_cu file
-!             LS convergence imposed from arm_cu file
+!             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
 !
       if (forcing_type .eq.0) THEN
@@ -269 +292 @@
       elseif (forcing_type .eq.5) THEN
        forcing_rico = .true.
+      elseif (forcing_type .eq.6) THEN
+       forcing_amma = .true.
       elseif (forcing_type .eq.40) THEN
        forcing_GCSSold = .true.
+      elseif (forcing_type .eq.59) THEN
+       forcing_sandu   = .true.
+      elseif (forcing_type .eq.60) THEN
+       forcing_astex   = .true.
       elseif (forcing_type .eq.61) THEN
        forcing_armcu = .true.
@@ -276 +305 @@
       else
        write (*,*) 'ERROR : unknown forcing_type ', forcing_type
-       stop 'Forcing_type should be 0,1,2,3 or 40'
+       stop 'Forcing_type should be 0,1,2,3,4,5,6 or 40,59,60,61'
       ENDIF
       print*,"forcing type=",forcing_type
@@ -286 +315 @@
 
         type_ts_forcing = 0
-        if (forcing_toga) type_ts_forcing = 1
+        if (forcing_toga .or. forcing_sandu .or. forcing_astex) 
+     :    type_ts_forcing = 1
 
 !---------------------------------------------------------------------
@@ -325 +355 @@
 c Special case for arm_cu which lasts less than one day : 53100s !! (MPL 20111026)
       IF(forcing_type .EQ. 61) fnday=53100./86400.
+c Special case for amma which lasts less than one day : 64800s !! (MPL 20120216)
+      IF(forcing_type .EQ. 6) fnday=64800./86400.
       annee_ref = anneeref
       mois = 1
@@ -334 +366 @@
       day_ini = day
       day_end = day_ini + nday
+
+      IF (forcing_type .eq.2) THEN
 ! Convert the initial date of Toga-Coare to Julian day
       call ymds2ju
      $ (year_ini_toga,mth_ini_toga,day_ini_toga,heure,day_ju_ini_toga)
 
+      ELSEIF (forcing_type .eq.4) THEN
 ! Convert the initial date of TWPICE to Julian day
       call ymds2ju
      $ (year_ini_twpi,mth_ini_twpi,day_ini_twpi,heure_ini_twpi
      $ ,day_ju_ini_twpi)
-
-! Convert the initial date of Arm_cu to Julian day
+      ELSEIF (forcing_type .eq.6) THEN
+! Convert the initial date of AMMA to Julian day
+      call ymds2ju
+     $ (year_ini_amma,mth_ini_amma,day_ini_amma,heure_ini_amma
+     $ ,day_ju_ini_amma)
+
+      ELSEIF (forcing_type .eq.59) THEN
+! Convert the initial date of Sandu case to Julian day
+      call ymds2ju
+     $   (year_ini_sandu,mth_ini_sandu,day_ini_sandu,
+     $    time_ini*3600.,day_ju_ini_sandu)
+
+      ELSEIF (forcing_type .eq.60) THEN
+! Convert the initial date of Astex case to Julian day
+      call ymds2ju
+     $   (year_ini_astex,mth_ini_astex,day_ini_astex,
+     $    time_ini*3600.,day_ju_ini_astex)
+
+      ELSEIF (forcing_type .eq.61) THEN
+
+! Convert the initial date of Arm_cu case to Julian day
       call ymds2ju
      $ (year_ini_armcu,mth_ini_armcu,day_ini_armcu,heure_ini_armcu
      $ ,day_ju_ini_armcu)
+      ENDIF
 
       daytime = day + time_ini/24. ! 1st day and initial time of the simulation 
@@ -418 +473 @@
 !!      preff= 1.01325e5
       preff = psurf
-      call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
+      IF (ok_old_disvert) THEN
+        call disvert0(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig)
+        print *,'On utilise disvert0'
+      ELSE
+        call disvert(pa,preff,ap,bp,dpres,presnivs,nivsigs,nivsig,
+     :                 scaleheight)
+        print *,'On utilise disvert'
+c       Nouvelle version disvert permettant d imposer ap,bp (modif L.Guez) MPL 18092012
+c       Dans ce cas, on lit ap,bp dans le fichier hybrid.txt
+      ENDIF
       sig_s=presnivs/preff
       plev =ap+bp*psurf
@@ -424 +488 @@
 ccc      zlay=-rd*300.*log(play/psurf)/rg ! moved after reading profiles
 
+      IF (forcing_type .eq. 59) THEN
+! pour forcing_sandu, on cherche l'indice le plus proche de 700hpa#3000m
       write(*,*) '***********************'
       do l = 1, llm
        write(*,*) 'l,play(l),presnivs(l): ',l,play(l),presnivs(l)
+       if (trouve_700 .and. play(l).le.70000) then
+         llm700=l
+         print *,'llm700,play=',llm700,play(l)/100.
+         trouve_700= .false.
+       endif
       enddo
       write(*,*) '***********************'
+      ENDIF
 
 c
@@ -460 +532 @@
 ! rday: defini dans suphel.F (86400.)
 ! day_ini: lu dans run.def (dayref)
-! rlat,rlon lus dans lmdz1d.def
+! rlat_rad,rlon-rad: transformes en radian de rlat,rlon lus dans lmdz1d.def (en degres)
 ! airefi,zcufi,zcvfi initialises au debut de ce programme
 ! rday,ra,rg,rd,rcpd declares dans YOMCST.h et calcules dans suphel.F
@@ -470 +542 @@
       zcufi=airefi
       zcvfi=airefi
+!
+      rlat_rad(:)=rlat(:)*rpi/180.
+      rlon_rad(:)=rlon(:)*rpi/180.
 
       call iniphysiq(ngrid,llm,rday,day_ini,timestep,
-     .     rlat,rlon,airefi,zcufi,zcvfi,ra,rg,rd,rcpd,1)
+     .     rlat_rad,rlon_rad,airefi,zcufi,zcvfi,ra,rg,rd,rcpd,(/1/))
       print*,'apres iniphysiq'
 
@@ -501 +576 @@
         agesno  = xagesno
         tsoil(:,:,:)=tsurf
+!------ AMMA 2e run avec modele sol et rayonnement actif (MPL 23052012)
+!       tsoil(1,1,1)=299.18
+!       tsoil(1,2,1)=300.08
+!       tsoil(1,3,1)=301.88
+!       tsoil(1,4,1)=305.48
+!       tsoil(1,5,1)=308.00
+!       tsoil(1,6,1)=308.00
+!       tsoil(1,7,1)=308.00
+!       tsoil(1,8,1)=308.00
+!       tsoil(1,9,1)=308.00
+!       tsoil(1,10,1)=308.00
+!       tsoil(1,11,1)=308.00
+!-----------------------------------------------------------------------
         call pbl_surface_init(qsol, fder, snsrf, qsurfsrf,
      &                                    evap, frugs, agesno, tsoil)
@@ -734 +822 @@
        endif
 
-       if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice) then
+       if (forcing_toga .or. forcing_GCSSold .or. forcing_twpice
+     :    .or.forcing_amma) then
          fcoriolis=0.0 ; ug=0. ; vg=0.
        endif
@@ -748 +837 @@
      : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
 
+!!!!!!!!!!!!!!!!!!!!!!!!
+! Geostrophic wind
+!!!!!!!!!!!!!!!!!!!!!!!!
+       sfdt = sin(0.5*fcoriolis*timestep)
+       cfdt = cos(0.5*fcoriolis*timestep)
+!
+        du_age(1:mxcalc)= -2.*sfdt/timestep*
+     :          (sfdt*(u(1:mxcalc)-ug(1:mxcalc)) -
+     :           cfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
+!!     : fcoriolis*(v(1:mxcalc)-vg(1:mxcalc))
+!
+       dv_age(1:mxcalc)= -2.*sfdt/timestep*
+     :          (cfdt*(u(1:mxcalc)-ug(1:mxcalc)) +
+     :           sfdt*(v(1:mxcalc)-vg(1:mxcalc))  )
+!!     : -fcoriolis*(u(1:mxcalc)-ug(1:mxcalc))
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!         call  writefield_phy('dv_age' ,dv_age,llm)
+!         call  writefield_phy('du_age' ,du_age,llm)
+!         call  writefield_phy('du_phys' ,du_phys,llm)
+!         call  writefield_phy('u_tend' ,u,llm)
+!         call  writefield_phy('u_g' ,ug,llm)
+!
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!! Increment state variables
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         u(1:mxcalc)=u(1:mxcalc) + timestep*(
      :              du_phys(1:mxcalc)
@@ -773 +888 @@
 
         teta=temp*(pzero/play)**rkappa
+!
+!---------------------------------------------------------------------
+!   Nudge soil temperature if requested
+!---------------------------------------------------------------------
+
+      IF (nudge_tsoil) THEN
+       ftsoil(1,isoil_nudge,:) = ftsoil(1,isoil_nudge,:) 
+     .  -timestep/tau_soil_nudge*(ftsoil(1,isoil_nudge,:)-Tsoil_nudge)
+      ENDIF
 
 !---------------------------------------------------------------------

LMDZ5/branches/testing/libf/phy1d/ocean_forced_mod.F90

@@ -31 +31 @@
     USE limit_read_mod
     USE mod_grid_phy_lmdz
-    INCLUDE "indicesol.h"
+    USE indice_sol_mod
+!    INCLUDE "indicesol.h"
     INCLUDE "YOMCST.h"
 
@@ -145 +146 @@
     USE limit_read_mod
     USE fonte_neige_mod,  ONLY : fonte_neige
-
-    INCLUDE "indicesol.h"
+    USE indice_sol_mod
+
+!    INCLUDE "indicesol.h"
     INCLUDE "dimsoil.h"
     INCLUDE "YOMCST.h"

LMDZ5/branches/testing/libf/phy1d/surf_land_bucket_mod.F90

@@ -26 +26 @@
     USE mod_grid_phy_lmdz
     USE mod_phys_lmdz_para
+    USE indice_sol_mod
 !****************************************************************************************
 ! Bucket calculations for surface. 
 !
     INCLUDE "clesphys.h"
-    INCLUDE "indicesol.h"
+!    INCLUDE "indicesol.h"
     INCLUDE "dimsoil.h"
     INCLUDE "YOMCST.h"