Index: dynamico_lmdz/simple_physics/phyparam/param/iniphyparam.F
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/param/iniphyparam.F	(revision 4214)
+++ dynamico_lmdz/simple_physics/phyparam/param/iniphyparam.F	(revision 4215)
@@ -12,5 +12,5 @@
       use comsaison
       USE geometry_mod, ONLY : longitude,latitude,cell_area
-      USE phys_const, ONLY : rad,g,r,cpp,rcp,dtphys,unjours,mugaz
+      USE phys_const, ONLY : planet_rad,g,r,cpp,rcp,dtphys,unjours,mugaz
       USE planet, ONLY : coefir, coefvis
       USE astronomy
@@ -95,5 +95,5 @@
       print*,'Avant les getpar '
       CALL getpar('unjours',86400.  ,unjours,'unjours')
-      CALL getpar('rad',6400000.    ,rad,'rad')
+      CALL getpar('planet_rad',6.4e6,planet_rad,'planet_rad')
       CALL getpar('g',9.8           ,g,'g')
       CALL getpar('cpp',1004.       ,cpp,'cpp')
@@ -129,5 +129,5 @@
 
       write(lunout,*) 'unjours=',unjours
-      write(lunout,*) 'rad=',rad
+      write(lunout,*) 'planet_rad=',planet_rad
       write(lunout,*) 'g=',g
       write(lunout,*) 'cpp=',cpp
@@ -197,5 +197,5 @@
       coslon(:)=cos(long(:))
 
-      prad=rad
+      prad=planet_rad
       pg=g
       r=8.134/(mugaz*0.001)
Index: dynamico_lmdz/simple_physics/phyparam/param/phyparam.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/param/phyparam.F90	(revision 4214)
+++ dynamico_lmdz/simple_physics/phyparam/param/phyparam.F90	(revision 4215)
@@ -1,8 +1,12 @@
 MODULE phyparam_mod
+  USE callkeys
+  USE comgeomfi
   IMPLICIT NONE
   PRIVATE
   SAVE
 
-  REAL, PARAMETER :: solarcst=1370., stephan=5.67e-08
+  REAL, PARAMETER :: pi=2*ASIN(1.), solarcst=1370., stephan=5.67e-08, &
+       ps_rad=1.e5, height_scale=10000., ref_temp=285., &
+       capcal_nosoil=1e5
     
   REAL, ALLOCATABLE :: tsurf(:),tsoil(:,:),rnatur(:), &
@@ -31,21 +35,11 @@
        &            pw, &
        &            pdu,pdv,pdt,pdq,pdpsrf)
-
-    USE callkeys
-    USE comsaison
-    USE dimphy
-    USE comgeomfi
-    USE phys_const
-    USE planet
-    USE astronomy
+    USE phys_const, ONLY : g, rcp, r, unjours
+    USE surface,    ONLY : soil
     USE turbulence, ONLY : vdif
     USE convection, ONLY : convadj
-    USE solar, ONLY : solang, zenang, mucorr
-    USE radiative_sw, ONLY : sw
-    USE radiative_lw, ONLY : lw
-    USE surface
 
     !=======================================================================
-    !   Organisation of the physical parametrisations of the LMD 
+    !   Top routine of the physical parametrisations of the LMD 
     !   20 parameters GCM for planetary atmospheres.
     !   It includes:
@@ -72,4 +66,5 @@
          pphi(ngrid,nlayer),    & ! Geopotential at the middle of the layers (m2s-2)
          pphis(ngrid),          & ! surface geopotential (unused)
+         presnivs(nlayer),      & 
          pu(ngrid,nlayer),      & ! u component of the wind (ms-1)
          pv(ngrid,nlayer),      & ! v component of the wind (ms-1)
@@ -95,10 +90,8 @@
     REAL zdhfr(ngrid,nlayer),zdtsrf(ngrid),zdtsrfr(ngrid)
     REAL zflubid(ngrid),zpmer(ngrid)
-    REAL zplanck(ngrid),zpopsk(ngrid,nlayer)
+    REAL zpopsk(ngrid,nlayer)
     REAL zdum1(ngrid,nlayer)
     REAL zdum2(ngrid,nlayer)
     REAL zdum3(ngrid,nlayer)
-    REAL ztim1,ztim2,ztim3
-    REAL zls,zinsol
     REAL zdtlw(ngrid,nlayer),zdtsw(ngrid,nlayer)
     REAL zfluxsw(ngrid),zfluxlw(ngrid)
@@ -109,11 +102,4 @@
     !   ----------------------
     
-    REAL zps_av
-    
-    INTEGER        longcles
-    PARAMETER    ( longcles = 20 )
-    REAL clesphy0( longcles      )
-    REAL presnivs(nlayer)
-        
     print*,'OK DANS PHYPARAM'
     print*,'nq ',nq
@@ -123,78 +109,4 @@
          &      size(pdq),size(lati),size(pq),size(factq)
     
-    nlevel=nlayer+1
-    igout=ngrid/2+1    
-    zday=rjourvrai+gmtime
-
-    pdq(:,:,:) = 0. ! we do not use tracers in this physics package
-
-    !-----------------------------------------------------------------------
-    !    1. Initialisations :
-    !    --------------------
-    
-    IF(firstcall) THEN
-       
-       print*,'AKk',ngrid,nsoilmx
-       allocate(tsurf(ngrid))
-       print*,'AKa'
-       allocate (tsoil(ngrid,nsoilmx))
-       print*,'AKb'
-       allocate (rnatur(ngrid))
-       print*,'AK2'
-       allocate(capcal(ngrid),fluxgrd(ngrid))
-       print*,'AK3'
-       allocate(dtrad(ngrid,nlayer),fluxrad(ngrid))
-       print*,'AK4'
-       allocate(q2(ngrid,nlayer+1),q2l(ngrid,nlayer+1))
-       print*,'AK5'
-       allocate(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid))
-       print*,'AK6'
-       
-       
-       PRINT*,'FIRSTCALL  '
-       
-       !         zday_last=rjourvrai
-       zday_last=zday-ptimestep/unjours
-       rnatur=1.
-       emissiv(:)=(1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter
-       inertie(:)=(1.-rnatur(:))*I_mer+rnatur(:)*I_ter
-       albedo(:)=(1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter
-       z0(:)=(1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter
-       q2=1.e-10
-       q2l=1.e-10
-       tsurf=300.
-       tsoil=300.
-       print*,tsoil(ngrid/2+1,nsoilmx/2+2)
-       print*,'TS ',tsurf(igout),tsoil(igout,5)
-       CALL iniorbit
-       
-       if (.not.callrad) then
-          do ig=1,ngrid
-             fluxrad(ig)=0.
-          enddo
-       endif
-       
-       !     print*,'OK PHYPARAM 1 '
-       IF(callsoil) THEN
-          CALL soil(ngrid,nsoilmx,firstcall,inertie, &
-          &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
-       ELSE
-          PRINT*,'WARNING!!! Thermal conduction in the soil turned off'
-          DO ig=1,ngrid
-             capcal(ig)=1.e5
-             fluxgrd(ig)=0.
-          ENDDO
-       ENDIF
-       !        CALL inifrict(ptimestep)
-       icount=0
-       
-       PRINT*,'FIRSTCALL B '
-       print*,'INIIO avant iophys_ini '
-       call iophys_ini('phys.nc    ',presnivs)
-       
-    ENDIF
-    icount=icount+1
-    
-    PRINT*,'FIRSTCALL AP '
     IF (ngrid.NE.ngridmax) THEN
        PRINT*,'STOP in inifis'
@@ -204,18 +116,50 @@
        STOP
     ENDIF
-    
-    DO l=1,nlayer
-       DO ig=1,ngrid
-          pdv(ig,l)=0.
-          pdu(ig,l)=0.
-          pdt(ig,l)=0.
-       ENDDO
-    ENDDO
-    
-    DO ig=1,ngrid
-       pdpsrf(ig)=0.
-       zflubid(ig)=0.
-       zdtsrf(ig)=0.
-    ENDDO
+
+    nlevel=nlayer+1
+    igout=ngrid/2+1    
+    zday=rjourvrai+gmtime
+
+    !-----------------------------------------------------------------------
+    !    0. Allocate and initialize at first call
+    !    --------------------
+    
+    IF(firstcall) THEN
+       !         zday_last=rjourvrai
+       zday_last=zday-ptimestep/unjours
+       CALL alloc_phyparam(ngrid, nlayer, igout)
+
+       !     print*,'OK PHYPARAM 1 '
+       IF(callsoil) THEN
+          CALL soil(ngrid,nsoilmx,firstcall,inertie, &
+               &          ptimestep,tsurf,tsoil,capcal,fluxgrd)
+          ! NB : this call to soil also performs some calculations, see surface.F90
+       ELSE
+          PRINT*,'WARNING!!! Thermal conduction in the soil turned off'
+          DO ig=1,ngrid
+             capcal(ig)  = capcal_nosoil
+             fluxgrd(ig) = 0.
+          ENDDO
+       ENDIF
+       !        CALL inifrict(ptimestep)
+       
+       PRINT*,'FIRSTCALL B '
+       print*,'INIIO avant iophys_ini '
+       call iophys_ini('phys.nc    ',presnivs)       
+    ENDIF
+
+    !-----------------------------------------------------------------------
+    !    1. Initialisations :
+    !    --------------------
+    
+    icount=icount+1
+    
+    pdq(:,:,:) = 0. ! we do not use tracers in this physics package    
+    pdv(:,:)  = 0.
+    pdu(:,:)  = 0.
+    pdt(:,:)  = 0.
+    pdpsrf(:) = 0.
+    zflubid(:)= 0.
+    zdtsrf(:) = 0.
     
     !-----------------------------------------------------------------------
@@ -243,9 +187,5 @@
     DO l=1,nlayer
        DO ig=1,ngrid
-          zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp
-       ENDDO
-    ENDDO
-    DO l=1,nlayer
-       DO ig=1,ngrid
+          zpopsk(ig,l)=(pplay(ig,l)/pplev(ig,1))**rcp ! surface pressure is used as reference pressure
           zh(ig,l)=pt(ig,l)/zpopsk(ig,l)
        ENDDO
@@ -256,112 +196,7 @@
     !    ---------------------------------------
     
-    !        print*,'callrad0'
-    IF(callrad) THEN
-       !     print*,'callrad'
-       
-       !   WARNING !!! on calcule le ray a chaque appel
-       !        IF( MOD(icount,iradia).EQ.0) THEN
-       
-       CALL solarlong(zday,zls)
-       CALL orbite(zls,dist_sol,declin)
-       !      declin=0.
-       !     print*,'ATTENTIOn : pdeclin = 0','  L_s=',zls
-       print*,'diurnal=',diurnal
-       IF(diurnal) THEN
-          if ( 1.eq.1 ) then
-             ztim1=SIN(declin)
-             ztim2=COS(declin)*COS(2.*pi*(zday-.5))
-             ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
-             CALL solang(ngrid,sinlon,coslon,sinlat,coslat, &
-             &         ztim1,ztim2,ztim3,                   &
-             &         mu0,fract)
-          else
-             zdtime=ptimestep*float(iradia)
-             CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract)
-             print*,'ZENANG '
-          endif
-          
-          IF(lverbose) THEN
-             PRINT*,'day, declin, sinlon,coslon,sinlat,coslat'
-             PRINT*,zday, declin,                       &
-             &            sinlon(igout),coslon(igout),  &
-             &            sinlat(igout),coslat(igout)
-          ENDIF
-       ELSE
-          print*,'declin,ngrid,rad',declin,ngrid,rad
-          CALL mucorr(ngrid,declin,lati,mu0,fract,10000.,rad)
-       ENDIF
-       
-       !    2.2 Calcul of the radiative tendencies and fluxes:
-       !    --------------------------------------------------
-       
-       !  2.1.2 levels
-       
-       zinsol=solarcst/(dist_sol*dist_sol)
-       zps_av=1.e5
-       if (firstcall) then
-          print*,'WARNING ps_rad impose'
-       endif
-       CALL sw(ngrid,nlayer,diurnal,coefvis,albedo, &
-       &              pplev,zps_av,                 &
-       &              mu0,fract,zinsol,             &
-       &              zfluxsw,zdtsw,                &
-       &              lverbose)
-       
-       CALL lw(ngrid,nlayer,coefir,emissiv, &
-       &             pplev,zps_av,tsurf,pt, &
-       &             zfluxlw,zdtlw,         &
-       &             lverbose)
-
-
-       !    2.4 total flux and tendencies:
-       !    ------------------------------
-       
-       !    2.4.1 fluxes
-       
-       DO ig=1,ngrid
-          fluxrad(ig)=emissiv(ig)*zfluxlw(ig) &
-               &         +zfluxsw(ig)*(1.-albedo(ig))
-          zplanck(ig)=tsurf(ig)*tsurf(ig)
-          zplanck(ig)=emissiv(ig)* &
-          &         stephan*zplanck(ig)*zplanck(ig)
-          fluxrad(ig)=fluxrad(ig)-zplanck(ig)
-       ENDDO
-       
-       !    2.4.2 temperature tendencies
-       
-       DO l=1,nlayer
-          DO ig=1,ngrid
-             dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
-          ENDDO
-       ENDDO
-       
-       !        ENDIF
-       
-       !    2.5 Transformation of the radiative tendencies:
-       !    -----------------------------------------------
-       
-       DO l=1,nlayer
-          DO ig=1,ngrid
-             pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
-          ENDDO
-       ENDDO
-       
-       IF(lverbose) THEN
-          PRINT*,'Diagnotique for the radaition'
-          PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'
-          PRINT*,albedo(igout),emissiv(igout),mu0(igout), &
-          &           fract(igout),                       &
-          &           fluxrad(igout),zplanck(igout)
-          PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'
-          PRINT*,'unjours',unjours
-          DO l=1,nlayer
-             WRITE(*,'(3f15.5,2e15.2)') pt(igout,l),   &
-             &         pplay(igout,l),pplev(igout,l),  &
-             &         zdtsw(igout,l),zdtlw(igout,l)
-          ENDDO
-       ENDIF 
-       
-    ENDIF
+    IF(callrad) CALL radiative_tendencies(ngrid, igout, nlayer, &
+         gmtime, ptimestep*float(iradia), zday, pplev, pplay, pt, &
+         pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, mu0)
 
     !-----------------------------------------------------------------------
@@ -444,5 +279,5 @@
 
     !-----------------------------------------------------------------------
-    !   On incremente les tendances physiques de la temperature du sol:
+    !   On ajoute les tendances physiques a la temperature du sol:
     !   ---------------------------------------------------------------
     
@@ -478,5 +313,5 @@
        
        do ig=1,ngridmax
-          zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(r*285.))
+          zpmer(ig)=pplev(ig,1)*exp(pphi(ig,1)/(r*ref_temp))
        enddo
        
@@ -523,6 +358,154 @@
   END SUBROUTINE phyparam
 
-
-  SUBROUTINE alloc_phyparam
+  SUBROUTINE radiative_tendencies(ngrid, igout, nlayer, &
+       gmtime, zdtime, zday, pplev, pplay, pt, &
+       pdt, zdtlw, zfluxlw, zdtsw, zfluxsw, mu0)
+    USE comsaison
+    USE planet
+    USE phys_const,   ONLY : planet_rad, unjours
+    USE astronomy,    ONLY : orbite, solarlong
+    USE solar,        ONLY : solang, zenang, mucorr
+    USE radiative_sw, ONLY : sw
+    USE radiative_lw, ONLY : lw
+
+    INTEGER, INTENT(IN) :: ngrid, igout, nlayer
+    REAL, INTENT(IN)    :: gmtime, zdtime, zday, &
+         &                 pplev(ngrid,nlayer+1), pplay(ngrid, nlayer), &
+         &                 pt(ngrid, nlayer+1)
+    REAL, INTENT(INOUT) :: pdt(ngrid,nlayer)
+    REAL, INTENT(OUT) :: zdtlw(ngrid,nlayer), zfluxlw(ngrid), &
+         zdtsw(ngrid,nlayer), zfluxsw(ngrid), mu0(ngrid)
+    REAL, DIMENSION(ngrid) :: fract
+    REAL :: zls, zinsol, tsurf2, ztim1,ztim2,ztim3
+    REAL :: zplanck(ngrid)
+    INTEGER :: ig, l
+
+    !    2.1 Insolation
+    !    --------------------------------------------------
+    
+    CALL solarlong(zday,zls)
+    CALL orbite(zls,dist_sol,declin)
+    
+    IF(diurnal) THEN
+       IF ( .TRUE. ) then
+          ztim1=SIN(declin)
+          ztim2=COS(declin)*COS(2.*pi*(zday-.5))
+          ztim3=-COS(declin)*SIN(2.*pi*(zday-.5))
+          CALL solang(ngrid,sinlon,coslon,sinlat,coslat, &
+               &         ztim1,ztim2,ztim3,                   &
+               &         mu0,fract)
+       ELSE
+          CALL zenang(ngrid,zls,gmtime,zdtime,lati,long,mu0,fract)
+          print*,'ZENANG '
+       ENDIF
+       
+       IF(lverbose) THEN
+          PRINT*,'day, declin, sinlon,coslon,sinlat,coslat'
+          PRINT*,zday, declin,                       &
+               &            sinlon(igout),coslon(igout),  &
+               &            sinlat(igout),coslat(igout)
+       ENDIF
+    ELSE
+       print*,'declin,ngrid,planet_rad',declin,ngrid,planet_rad
+       CALL mucorr(ngrid,declin,lati,mu0,fract,height_scale,planet_rad)
+    ENDIF
+    
+    zinsol=solarcst/(dist_sol*dist_sol)
+    
+    !    2.2 Radiative tendencies and fluxes:
+    !    --------------------------------------------------
+    
+    CALL sw(ngrid,nlayer,diurnal,coefvis,albedo, &
+         &              pplev,ps_rad,                 &
+         &              mu0,fract,zinsol,             &
+         &              zfluxsw,zdtsw,                &
+         &              lverbose)
+    
+    CALL lw(ngrid,nlayer,coefir,emissiv, &
+         &             pplev,ps_rad,tsurf,pt, &
+         &             zfluxlw,zdtlw,         &
+         &             lverbose)
+    
+    !    2.4 surface fluxes
+    !    ------------------------------
+    
+    DO ig=1,ngrid
+       fluxrad(ig)=emissiv(ig)*zfluxlw(ig) &
+            &         +zfluxsw(ig)*(1.-albedo(ig))
+       tsurf2 = tsurf(ig)*tsurf(ig)
+       zplanck(ig)=emissiv(ig)*stephan*tsurf2*tsurf2
+       fluxrad(ig)=fluxrad(ig)-zplanck(ig)
+    ENDDO
+    
+    !    2.5 Temperature tendencies
+    !    --------------------------
+    
+    DO l=1,nlayer
+       DO ig=1,ngrid
+          dtrad(ig,l)=zdtsw(ig,l)+zdtlw(ig,l)
+          pdt(ig,l)=pdt(ig,l)+dtrad(ig,l)
+       ENDDO
+    ENDDO
+    
+    IF(lverbose) THEN
+       PRINT*,'Diagnostics for radiation'
+       PRINT*,'albedo, emissiv, mu0,fract,Frad,Planck'
+       PRINT*,albedo(igout),emissiv(igout),mu0(igout), &
+            &           fract(igout),                       &
+            &           fluxrad(igout),zplanck(igout)
+       PRINT*,'Tlay Play Plev dT/dt SW dT/dt LW (K/day)'
+       PRINT*,'unjours',unjours
+       DO l=1,nlayer
+          WRITE(*,'(3f15.5,2e15.2)') pt(igout,l),   &
+               &         pplay(igout,l),pplev(igout,l),  &
+               &         zdtsw(igout,l),zdtlw(igout,l)
+       ENDDO
+    ENDIF
+    
+  END SUBROUTINE radiative_tendencies
+
+  SUBROUTINE alloc_phyparam(ngrid, nlayer, igout)
+    USE surface
+    USE astronomy, ONLY : iniorbit
+    INTEGER, INTENT(IN) :: ngrid, nlayer, igout
+    LOGICAL, PARAMETER :: firstcall=.TRUE.
+
+    print*,'AKk',ngrid,nsoilmx
+    allocate(tsurf(ngrid))
+    print*,'AKa'
+    allocate (tsoil(ngrid,nsoilmx))
+    print*,'AKb'
+    allocate (rnatur(ngrid))
+    print*,'AK2'
+    allocate(capcal(ngrid),fluxgrd(ngrid))
+    print*,'AK3'
+    allocate(dtrad(ngrid,nlayer),fluxrad(ngrid))
+    print*,'AK4'
+    allocate(q2(ngrid,nlayer+1),q2l(ngrid,nlayer+1))
+    print*,'AK5'
+    allocate(albedo(ngrid),emissiv(ngrid),z0(ngrid),inertie(ngrid))
+    print*,'AK6'
+    
+    
+    PRINT*,'FIRSTCALL  '
+    
+    rnatur=1.
+    emissiv(:)=(1.-rnatur(:))*emi_mer+rnatur(:)*emi_ter
+    inertie(:)=(1.-rnatur(:))*I_mer+rnatur(:)*I_ter
+    albedo(:)=(1.-rnatur(:))*alb_mer+rnatur(:)*alb_ter
+    z0(:)=(1.-rnatur(:))*Cd_mer+rnatur(:)*Cd_ter
+    q2=1.e-10
+    q2l=1.e-10
+    tsurf=300.
+    tsoil=300.
+    print*,tsoil(igout,nsoilmx/2+2)
+    print*,'TS ',tsurf(igout),tsoil(igout,5)
+    CALL iniorbit
+    
+    if (.not.callrad) fluxrad(:)=0.
+    
+    icount=0
+    
   END SUBROUTINE alloc_phyparam
+
 END MODULE phyparam_mod
Index: dynamico_lmdz/simple_physics/phyparam/physics/phys_const.F90
===================================================================
--- dynamico_lmdz/simple_physics/phyparam/physics/phys_const.F90	(revision 4214)
+++ dynamico_lmdz/simple_physics/phyparam/physics/phys_const.F90	(revision 4215)
@@ -1,5 +1,5 @@
 MODULE phys_const
   
-  REAL rad,g,r,cpp,rcp,dtphys,unjours,mugaz,omeg
+  REAL planet_rad,g,r,cpp,rcp,dtphys,unjours,mugaz,omeg
 
 END MODULE phys_const
