Index: trunk/LMDZ.VENUS/libf/phyvenus/physiq.F
===================================================================
--- trunk/LMDZ.VENUS/libf/phyvenus/physiq.F	(revision 1306)
+++ trunk/LMDZ.VENUS/libf/phyvenus/physiq.F	(revision 1310)
@@ -56,5 +56,4 @@
       USE ioipsl
 !      USE histcom ! not needed; histcom is included in ioipsl
-      USE chemparam_mod
       USE infotrac
       USE control_mod
@@ -66,4 +65,8 @@
       USE iophy
       use cpdet_mod, only: cpdet, t2tpot
+      USE chemparam_mod
+      USE conc
+      USE compo_hedin83_mod2
+      use moyzon_mod, only: tmoy
       use ieee_arithmetic
       IMPLICIT none
@@ -86,4 +89,6 @@
 #include "logic.h"
 #include "tabcontrol.h"
+#include "nirdata.h"
+#include "hedin.h"
 c======================================================================
       LOGICAL ok_journe ! sortir le fichier journalier
@@ -223,4 +228,16 @@
       EXTERNAL mucorr
       EXTERNAL phytrac
+      EXTERNAL nirco2abs
+      EXTERNAL nir_leedat
+      EXTERNAL nltecool
+      EXTERNAL nlte_tcool
+      EXTERNAL nlte_setup
+      EXTERNAL blendrad
+      EXTERNAL nlthermeq
+      EXTERNAL euvheat
+      EXTERNAL param_read
+      EXTERNAL param_read_e107
+      EXTERNAL conduction
+      EXTERNAL molvis
 c
 c Variables locales
@@ -247,4 +264,9 @@
       REAL zphi(klon,klev)
       REAL zzlev(klon,klev+1),zzlay(klon,klev),z1,z2
+      real tsurf(klon)
+
+c va avec nlte_tcool
+      INTEGER ierr_nlte
+      REAL    varerr
 
 c Variables du changement
@@ -268,4 +290,14 @@
       REAL d_u_hin(klon,klev), d_v_hin(klon,klev)
       REAL d_t_hin(klon,klev)
+
+c Tendencies due to radiative scheme   [K/s]
+c     d_t_rad,dtsw,dtlw,d_t_nirco2,d_t_nlte,d_t_euv
+c are not computed at each physical timestep
+c therefore, they are defined and saved in phys_state_var_mod
+
+c Tendencies due to molecular viscosity and conduction
+      real d_t_conduc(klon,klev)     ! [K/s]
+      real d_u_molvis(klon,klev)     ! (m/s) /s
+      real d_v_molvis(klon,klev)     ! (m/s) /s
 
 c
@@ -386,4 +418,8 @@
          call phys_state_var_init
 c
+c Initialising Hedin model for upper atm 
+c   (to be revised when coupled to chemistry) :
+         call conc_init
+c
 c Initialiser les compteurs:
 c
@@ -458,6 +494,36 @@
 C source dans couche limite
          source = 0.0 ! pas de source, pour l'instant
-C
 c---------
+
+c---------
+c INITIALIZE THERMOSPHERIC PARAMETERS
+c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+         if (callthermos) then
+            if(solvarmod.eq.0) call param_read
+            if(solvarmod.eq.1) call param_read_e107 
+         endif
+
+c Initialisation (recomputed in concentration2)
+       do ig=1,klon
+         do j=1,klev
+            rnew(ig,j)=R
+            cpnew(ig,j)=cpdet(tmoy(j))
+c            print*, ' physique  l503'
+c            print*,  j, cpdet(tmoy(j))
+             mmean(ig,j)=RMD
+            akknew(ig,j)=1.e-4
+          enddo
+c        stop
+
+        enddo  
+     
+      IF(callthermos.or.callnlte.or.callnirco2) THEN  
+         call compo_hedin83_init2
+      ENDIF
+      if (callnlte) call nlte_setup
+      if(callnirco2.and.(nircorr.eq.1)) call nir_leedat         
+c---------
+      
 c
 c Verifications:
@@ -684,9 +750,9 @@
       DO k=1,klev
          DO i=1,klon
-            zzlay(i,k)=zphi(i,k)/RG
+            zzlay(i,k)=zphi(i,k)/RG        ! [m]
          ENDDO
       ENDDO
       DO i=1,klon
-         zzlev(i,1)=pphis(i)/RG
+         zzlev(i,1)=pphis(i)/RG            ! [m]
       ENDDO
       DO k=2,klev
@@ -1042,38 +1108,171 @@
 c====================================================================
 
+c-----------------------------------------------------------------------
+c    . Compute radiative tendencies :
+c------------------------------------
+c====================================================================
       IF (MOD(itaprad,radpas).EQ.0) THEN
-c             print*,'RAYONNEMENT ', 
-c    $             ' (itaprad=',itaprad,'/radpas=',radpas,')'
+c====================================================================
 
       dtimerad = dtime*REAL(radpas)  ! pas de temps du rayonnement (s)
 c     PRINT*,'dtimerad,dtime,radpas',dtimerad,dtime,radpas
             
+c Calcul pour Cp rnew et mmean avec traceurs (Cp independant de T !! )
+       IF(callnlte.or.callthermos) THEN                                 
+         call compo_hedin83_mod(pplay,rmu0,   
+     &	         co2vmr_gcm,covmr_gcm,ovmr_gcm,n2vmr_gcm,nvmr_gcm)
+       ENDIF	
+
+      if(callthermos) then
+         call concentrations2(pplay,t_seri,d_t,co2vmr_gcm, n2vmr_gcm,
+     &		covmr_gcm, ovmr_gcm,nvmr_gcm,pdtphys)
+      endif
+
+
+c          NLTE cooling from CO2 emission
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        IF(callnlte) THEN                                 
+            if(nltemodel.eq.0.or.nltemodel.eq.1) then
+                CALL nltecool(klon, klev, pplay*9.869e-6, t_seri,
+     $                    co2vmr_gcm,n2vmr_gcm, covmr_gcm, ovmr_gcm, 
+     $                    d_t_nlte)
+            else if(nltemodel.eq.2) then                                
+                CALL nlte_tcool(klon,klev,pplay*9.869e-6,              
+     $               t_seri,zzlay,co2vmr_gcm, n2vmr_gcm, covmr_gcm, 
+     $               ovmr_gcm,d_t_nlte,ierr_nlte,varerr )
+                  if(ierr_nlte.gt.0) then
+                     write(*,*)
+     $                'WARNING: nlte_tcool output with error message',
+     $                'ierr_nlte=',ierr_nlte,'varerr=',varerr
+                     write(*,*)'I will continue anyway'
+                  endif
+
+             endif
+             
+        ELSE
+ 
+          d_t_nlte(:,:)=0.
+
+        ENDIF	
+
+c      Find number of layers for LTE radiation calculations
+
+      IF(callnlte .or. callnirco2) 
+     $        CALL nlthermeq(klon, klev, paprs, pplay) 
+
+
+c          LTE radiative transfert / solar / IR matrix
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
       CALL radlwsw
      e            (dist, rmu0, fract, zzlev,
-     e             paprs, pplay,ftsol, t_seri,
-     s             heat,cool,radsol,
-     s             topsw,toplw,solsw,sollw,
-     s             sollwdown,
-     s             lwnet, swnet)
-
-      itaprad = 0
+     e             paprs, pplay,ftsol, t_seri)
+
+
+c          CO2 near infrared absorption
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        d_t_nirco2(:,:)=0.
+        if (callnirco2) then
+           call nirco2abs (klon, klev, pplay, dist, nqmax, qx,
+     .                 rmu0, fract, d_t_nirco2,
+     .                 co2vmr_gcm, ovmr_gcm)
+        endif
+
+
+c          Net atmospheric radiative heating rate (K.s-1)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        IF(callnlte.or.callnirco2) THEN 
+           CALL blendrad(klon, klev, pplay,heat,
+     &          cool, d_t_nirco2,d_t_nlte, dtsw, dtlw)
+        ELSE
+	   dtsw(:,:)=heat(:,:)
+	   dtlw(:,:)=-1*cool(:,:)
+	ENDIF
+
+         DO k=1,klev
+            DO i=1,klon
+               d_t_rad(i,k) = dtsw(i,k) + dtlw(i,k)   ! K/s
+            ENDDO
+         ENDDO
+
+
+cc---------------------------------------------
+
+c          EUV heating rate (K.s-1)
+c          ~~~~~~~~~~~~~~~~~~~~~~~~
+
+        d_t_euv(:,:)=0.
+
+        IF (callthermos) THEN
+
+c        call thermosphere(zplev,zplay,dist_sol,
+c     $     mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay,
+c     &     pt,pq,pu,pv,pdt,pdq,
+c     $     zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff)
+
+           call euvheat(klon, klev, t_seri,paprs,pplay,zzlay,
+     $          rmu0,pdtphys,gmtime,rjourvrai,
+C     $			pq,pdq,zdteuv)
+     $          co2vmr_gcm, n2vmr_gcm, covmr_gcm, 
+     $          ovmr_gcm,nvmr_gcm,d_t_euv )
+
+           DO k=1,klev
+              DO ig=1,klon
+                 d_t_rad(ig,k)=d_t_rad(ig,k)+d_t_euv(ig,k)
+                
+              ENDDO
+           ENDDO
+
+        ENDIF  ! callthermos
+
+c====================================================================
+        itaprad = 0
+      ENDIF    ! radpas
+c====================================================================
+c
+c Ajouter la tendance des rayonnements (tous les pas)
+c
       DO k = 1, klev
       DO i = 1, klon
-         dtrad(i,k) = heat(i,k)-cool(i,k)  !K/s
-      ENDDO
-      ENDDO
-
-      ENDIF
-      itaprad = itaprad + 1
-c====================================================================
-c
-c Ajouter la tendance des rayonnements (tous les pas)
-c
-      DO k = 1, klev
-      DO i = 1, klon
-         t_seri(i,k) = t_seri(i,k) + dtrad(i,k) * dtime
-      ENDDO
-      ENDDO
-
+         t_seri(i,k) = t_seri(i,k) + d_t_rad(i,k) * dtime
+      ENDDO
+      ENDDO
+
+! CONDUCTION  and  MOLECULAR VISCOSITY
+c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+        d_t_conduc(:,:)=0.
+        d_u_molvis(:,:)=0.
+        d_v_molvis(:,:)=0.
+
+        IF (callthermos) THEN
+
+           tsurf(:)=t_seri(:,1)
+           call conduction(klon, klev,pdtphys,
+     $            pplay,paprs,t_seri,
+     $            tsurf,zzlev,zzlay,d_t_conduc)
+
+            call molvis(klon, klev,pdtphys,
+     $            pplay,paprs,t_seri,
+     $            u,tsurf,zzlev,zzlay,d_u_molvis)
+
+            call molvis(klon, klev, pdtphys,
+     $            pplay,paprs,t_seri,
+     $            v,tsurf,zzlev,zzlay,d_u_molvis)
+
+            DO k=1,klev
+               DO ig=1,klon
+                  t_seri(ig,k)= t_seri(ig,k)+ d_t_conduc(ig,k)*dtime ! [K]
+                  u_seri(ig,k)= u_seri(ig,k)+ d_u_molvis(ig,k)*dtime ! m/s
+                  v_seri(ig,k)= v_seri(ig,k)+ d_v_molvis(ig,k)*dtime ! m/s
+               ENDDO
+            ENDDO
+
+         ENDIF  ! callthermos
+
+c====================================================================
       ! tests: output tendencies
 !      call writefield_phy('physiq_dtrad',dtrad,klev)
