Index: trunk/libf/dyn3dpar/calfis_p.F
===================================================================
--- trunk/libf/dyn3dpar/calfis_p.F	(revision 37)
+++ trunk/libf/dyn3dpar/calfis_p.F	(revision 52)
@@ -288,5 +288,5 @@
       ALLOCATE(flxwfi(klon,llm))
 ! ADAPTATION GCM POUR CP(T)
-      ALLOCATE(zteta(klon,llm)
+      ALLOCATE(zteta(klon,llm))
       ALLOCATE(zpk(klon,llm))
 c$OMP END MASTER
Index: trunk/libf/dyn3dpar/comconst.h
===================================================================
--- trunk/libf/dyn3dpar/comconst.h	(revision 37)
+++ trunk/libf/dyn3dpar/comconst.h	(revision 52)
@@ -11,5 +11,5 @@
      &                   ,dissip_factz,dissip_deltaz,dissip_zref        &
      &                   ,tau_top_bound,                                &
-     & daylen,year_day,molmass
+     & daylen,year_day,molmass, ihf
       COMMON/cpdetvenus/nu_venus,t0_venus
 
@@ -36,4 +36,5 @@
 
       REAL nu_venus,t0_venus ! coeffs needed for Cp(T), Venus atmosphere
+      REAL ihf  ! (W/m2) intrinsic heat flux for giant planets
 
 
Index: trunk/libf/dyn3dpar/conf_planete.F90
===================================================================
--- trunk/libf/dyn3dpar/conf_planete.F90	(revision 37)
+++ trunk/libf/dyn3dpar/conf_planete.F90	(revision 52)
@@ -64,3 +64,11 @@
 CALL getin('omeg',omeg)
 
+! Intrinsic heat flux [default is none]
+! Aymeric -- for giant planets
+! [matters only if planet_type="giant"]
+ihf = 0.
+CALL getin('ihf',ihf)
+
+
+
 END SUBROUTINE conf_planete
Index: trunk/libf/dyn3dpar/gcm.F
===================================================================
--- trunk/libf/dyn3dpar/gcm.F	(revision 37)
+++ trunk/libf/dyn3dpar/gcm.F	(revision 52)
@@ -59,4 +59,5 @@
 c   Declarations:
 c   -------------
+
 #include "dimensions.h"
 #include "paramet.h"
@@ -67,4 +68,5 @@
 #include "logic.h"
 #include "temps.h"
+!!!!!!!!!!!#include "control.h"
 #include "ener.h"
 #include "description.h"
@@ -73,10 +75,8 @@
 #include "iniprint.h"
 #include "tracstoke.h"
-
 #ifdef INCA
 ! Only INCA needs these informations (from the Earth's physics)
 #include "indicesol.h"
 #endif
-
       INTEGER         longcles
       PARAMETER     ( longcles = 20 )
@@ -93,5 +93,5 @@
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
       REAL teta(ip1jmp1,llm)                 ! temperature potentielle 
-      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: q ! champs advectes
+      REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
       REAL ps(ip1jmp1)                       ! pression  au sol
 c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
@@ -219,8 +219,9 @@
 #endif
 !      endif ! of if (planet_type.eq."earth")
-
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+c
 c Initialisations pour Cp(T) Venus
       call ini_cpdet
-
+c
 c-----------------------------------------------------------------------
 c   Choix du calendrier
@@ -244,4 +245,5 @@
       endif
 #endif
+c-----------------------------------------------------------------------
 
       IF (config_inca /= 'none') THEN
@@ -289,5 +291,5 @@
      &              teta,q,masse,ps,phis, time_0)
         endif ! of if (planet_type.eq."mars")
-
+        
 c       write(73,*) 'ucov',ucov
 c       write(74,*) 'vcov',vcov
@@ -304,4 +306,5 @@
          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
       endif
+
 
 c-----------------------------------------------------------------------
@@ -350,4 +353,5 @@
         write(lunout,*)' Pas de remise a zero'
       ENDIF
+
 c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
 c        write(lunout,*)
@@ -432,5 +436,6 @@
 c   Initialisation de la physique :
 c   -------------------------------
-      IF (call_iniphys.and.iflag_phys.eq.1) THEN
+
+      IF (call_iniphys.and.(iflag_phys.eq.1)) THEN
          latfi(1)=rlatu(1)
          lonfi(1)=0.
@@ -450,7 +455,7 @@
          zcvfi(ngridmx) = cv(ip1jm-iim)
          CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
-
          WRITE(lunout,*)
      .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
+
 ! Initialisation de la physique: pose probleme quand on tourne 
 ! SANS physique, car iniphysiq.F est dans le repertoire phy[]... 
@@ -462,5 +467,5 @@
          call_iniphys=.false.
       ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
-
+!#endif
 
 c-----------------------------------------------------------------------
@@ -508,10 +513,9 @@
       IF (mpi_rank==0) then
         if (ok_dyn_ins) then
-          ! initialize output file for instantaneous outputs
-          ! t_ops = iecri * daysec ! do operations every t_ops
-          t_ops =((1.0*iecri)/day_step) * daysec  
-          t_wrt = daysec ! iecri * daysec ! write output every t_wrt
-          t_wrt = daysec ! iecri * daysec ! write output every t_wrt
-          CALL inithist(day_ref,annee_ref,time_step,
+        ! initialize output file for instantaneous outputs
+        ! t_ops = iecri * daysec ! do operations every t_ops
+        t_ops =((1.0*iecri)/day_step) * daysec  
+        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
+        CALL inithist(day_ref,annee_ref,time_step,
      &                  t_ops,t_wrt)
         endif
Index: trunk/libf/dyn3dpar/iniacademic.F
===================================================================
--- trunk/libf/dyn3dpar/iniacademic.F	(revision 37)
+++ trunk/libf/dyn3dpar/iniacademic.F	(revision 52)
@@ -96,4 +96,7 @@
 ! --------------------------------------
 c
+
+        print *, 'This is iniacademic' 
+
         ! initialize planet radius, rotation rate,...
         call conf_planete
@@ -155,4 +158,6 @@
           teta0=315.     ! mean Teta (S.H. 315K)
           CALL getin('teta0',teta0)
+          print *, 'iniacademic - teta0 ', teta0
+          print *, 'iniacademic - rad ', rad
           ttp=200.       ! Tropopause temperature (S.H. 200K)
           CALL getin('ttp',ttp)
@@ -200,4 +205,14 @@
              tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
      &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
+             !! Aymeric -- tests particuliers
+             if (planet_type=="giant") then
+             tetajl(j,l)=teta0+(delt_y*
+     &          ((sin(rlatu(j)*3.14159*eps+0.0001))**2)
+     &          / ((rlatu(j)*3.14159*eps+0.0001)**2))
+     &          -delt_z*log(zsig)
+!!!             ddsin=sin(2.5*3.14159*rlatu(j))
+!!!             tetajl(j,l)=teta0-delt_y*ddsin*ddsin
+!!!!     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
+             endif
              ! Profil stratospherique isotherme (+vortex)
              w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
@@ -217,4 +232,5 @@
             enddo
           enddo
+          PRINT *, 'iniacademic - check',tetajl(:,int(llm/2)),rlatu(:)
 
 
Index: trunk/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- trunk/libf/dyn3dpar/leapfrog_p.F	(revision 37)
+++ trunk/libf/dyn3dpar/leapfrog_p.F	(revision 52)
@@ -117,4 +117,8 @@
       REAL,DIMENSION(:,:,:),ALLOCATABLE,SAVE :: dqfi
 
+       !! Aymeric -- cp(T) comme dans leapfrog.F, SAVE OK ???
+      REAL,SAVE :: duspg(ip1jmp1,llm) ! for bilan_dyn
+
+
 c   variables pour le fichier histoire
       REAL dtav      ! intervalle de temps elementaire
@@ -177,5 +181,10 @@
 
       logical , parameter :: flag_verif = .false.
-      
+
+      ! for CP(T)  -- Aymeric
+      real :: dtec
+      real,external :: cpdet
+      real,save :: ztetaec(ip1jmp1,llm)  !!SAVE ???
+
 c declaration liees au parallelisme
       INTEGER :: ierr
@@ -560,5 +569,5 @@
 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       do l=1,llm
-        tsurpk(ijb:ije,llm)=cpp*temp(ijb:ije,llm)/pk(ijb:ije,llm)
+        tsurpk(ijb:ije,l)=cpp*temp(ijb:ije,l)/pk(ijb:ije,l)
       enddo
 !$OMP END DO
@@ -999,4 +1008,16 @@
        enddo ! of do l=1,llm
 !$OMP END DO
+
+       if (planet_type.eq."giant") then
+          ! Intrinsic heat flux
+          ! Aymeric -- for giant planets
+          if (ihf .gt. 1.e-6) then
+          !print *, '**** INTRINSIC HEAT FLUX ****', ihf 
+          teta(ijb:ije,1) = teta(ijb:ije,1)
+     &        + dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
+          !print *, '**** d teta ' 
+          !print *, dtvr * aire(ijb:ije) * ihf / cpp / masse(ijb:ije,1)
+          endif
+       endif
 
        call Register_Hallo(ucov,ip1jmp1,llm,0,1,1,0,Request_Physic)
@@ -1463,5 +1484,5 @@
 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
       do l=1,llm
-        tsurpk(ijb:ije,llm)=cpp*temp(ijb:ije,llm)/pk(ijb:ije,llm)
+        tsurpk(ijb:ije,l)=cpp*temp(ijb:ije,l)/pk(ijb:ije,l)
       enddo
 !$OMP END DO
@@ -1685,7 +1706,7 @@
                 ije=ij_end
 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
-                do l=1,llm
-                  tsurpk(ijb:ije,llm)=cpp*temp(ijb:ije,llm)/
-     &                                             pk(ijb:ije,llm)
+                do l=1,llm      
+                  tsurpk(ijb:ije,l)=cpp*temp(ijb:ije,l)/
+     &                                             pk(ijb:ije,l)
                 enddo
 !$OMP END DO
