Index: /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/academic.h
===================================================================
--- /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/academic.h	(revision 1436)
+++ /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/academic.h	(revision 1437)
@@ -1,5 +1,9 @@
 !
-! $Header$
+! $Id$
 !
-      real tetarappel(ip1jmp1,llm),taurappel
-      common/academic/tetarappel,taurappel
+      common/academic/tetarappel,knewt_t,kfrict,knewt_g,clat4
+      real :: tetarappel(ip1jmp1,llm)
+      real :: knewt_t(llm)
+      real :: kfrict(llm)
+      real :: knewt_g
+      real :: clat4(ip1jmp1)
Index: /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/comconst.h
===================================================================
--- /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/comconst.h	(revision 1436)
+++ /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/comconst.h	(revision 1437)
@@ -5,18 +5,33 @@
 ! INCLUDE comconst.h
 
-      COMMON/comconst/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,           &
-     & dtvr,daysec,                                                     &
+      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
+     &                 iflag_top_bound
+      COMMON/comconstr/dtvr,daysec,                                     &
      & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg              &
      &                   ,dissip_factz,dissip_deltaz,dissip_zref        &
-     &                   ,iflag_top_bound,tau_top_bound
+     &                   ,tau_top_bound,                                &
+     & daylen,year_day,molmass
 
 
       INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
-      REAL dtvr,daysec
-      REAL pi,dtphys,dtdiss,rad,r,cpp,kappa
-      REAL cotot,unsim,g,omeg
+      REAL dtvr ! dynamical time step (in s)
+      REAL daysec !length (in s) of a standard day
+      REAL pi    ! something like 3.14159....
+      REAL dtphys ! (s) time step for the physics
+      REAL dtdiss ! (s) time step for the dissipation
+      REAL rad ! (m) radius of the planet
+      REAL r ! Gas constant R=8.31 J.K-1.mol-1
+      REAL cpp   ! Cp
+      REAL kappa ! kappa=R/Cp 
+      REAL cotot
+      REAL unsim ! = 1./iim
+      REAL g ! (m/s2) gravity
+      REAL omeg ! (rad/s) rotation rate of the planet
       REAL dissip_factz,dissip_deltaz,dissip_zref
       INTEGER iflag_top_bound
       REAL tau_top_bound
+      REAL daylen ! length of solar day, in 'standard' day length
+      REAL year_day ! Number of standard days in a year
+      REAL molmass ! (g/mol) molar mass of the atmosphere
 
 
Index: /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/conf_planete.F90
===================================================================
--- /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/conf_planete.F90	(revision 1437)
+++ /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/conf_planete.F90	(revision 1437)
@@ -0,0 +1,66 @@
+!
+! $Id$
+!
+SUBROUTINE conf_planete
+!
+#ifdef CPP_IOIPSL
+USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+USE ioipsl_getincom
+#endif
+IMPLICIT NONE
+!
+!
+!   Declarations :
+!   --------------
+#include "dimensions.h"
+#include "comconst.h"
+#include "comvert.h"
+!
+!   local:
+!   ------
+
+! ---------------------------------------------
+! Initialisations de constantes de la dynamique
+! ---------------------------------------------
+! Pi
+pi=2.*asin(1.)
+
+!Reference surface pressure (Pa)
+preff=101325.
+CALL getin('preff', preff)
+! Reference pressure at which hybrid coord. become purely pressure
+! pa=50000.
+pa=preff/2.
+CALL getin('pa', pa)
+! Gravity
+g=9.80665
+CALL getin('g',g)
+! Molar mass of the atmosphere
+molmass = 28.9644
+CALL getin('molmass',molmass)
+! kappa=R/Cp et Cp      
+kappa = 2./7.
+CALL getin('kappa',kappa)
+cpp=8.3145/molmass/kappa*1000.
+CALL getin('cpp',cpp)
+! Radius of the planet
+rad = 6371229. 
+CALL getin('radius',rad)
+! Length of a standard day (s)
+daysec=86400.
+CALL getin('daysec',daysec)
+! Rotation rate of the planet:
+! Length of a solar day, in standard days
+daylen = 1.
+CALL getin('daylen',daylen)
+! Number of days (standard) per year:
+year_day = 365.25
+CALL getin('year_day',year_day)
+! Omega
+! omeg=2.*pi/86400.
+omeg=2.*pi/daysec*(1./daylen+1./year_day)
+CALL getin('omeg',omeg)
+
+END SUBROUTINE conf_planete
Index: /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/friction.F
===================================================================
--- /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/friction.F	(revision 1436)
+++ /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/friction.F	(revision 1437)
@@ -6,18 +6,24 @@
 
       USE control_mod
- 
+#ifdef CPP_IOIPSL
+      USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      USE ioipsl_getincom
+#endif
+      
       IMPLICIT NONE
 
-c=======================================================================
-c
-c
-c   Objet:
-c   ------
-c
-c  ***********
-c    Friction
-c  ***********
-c
-c=======================================================================
+!=======================================================================
+!
+!   Friction for the Newtonian case:
+!   --------------------------------
+!    2 possibilities (depending on flag 'friction_type'
+!     friction_type=0 : A friction that is only applied to the lowermost
+!                       atmospheric layer
+!     friction_type=1 : Friction applied on all atmospheric layer (but
+!       (default)       with stronger magnitude near the surface; see
+!                       iniacademic.F)
+!=======================================================================
 
 #include "dimensions.h"
@@ -25,15 +31,36 @@
 #include "comgeom2.h"
 #include "comconst.h"
+#include "iniprint.h"
+#include "academic.h"
 
-      REAL pdt
+! arguments:
+      REAL,INTENT(out) :: ucov( iip1,jjp1,llm )
+      REAL,INTENT(out) :: vcov( iip1,jjm,llm )
+      REAL,INTENT(in) :: pdt ! time step
+
+! local variables:
+
       REAL modv(iip1,jjp1),zco,zsi
       REAL vpn,vps,upoln,upols,vpols,vpoln
       REAL u2(iip1,jjp1),v2(iip1,jjm)
-      REAL ucov( iip1,jjp1,llm ),vcov( iip1,jjm,llm )
-      INTEGER  i,j
-      REAL cfric
-      parameter (cfric=1.e-5)
+      INTEGER  i,j,l
+      REAL,PARAMETER :: cfric=1.e-5
+      LOGICAL,SAVE :: firstcall=.true.
+      INTEGER,SAVE :: friction_type=1
+      CHARACTER(len=20) :: modname="friction"
+      CHARACTER(len=80) :: abort_message
+      
+      IF (firstcall) THEN
+        ! set friction type
+        call getin("friction_type",friction_type)
+        if ((friction_type.lt.0).or.(friction_type.gt.1)) then
+          abort_message="wrong friction type"
+          write(lunout,*)'Friction: wrong friction type',friction_type
+          call abort_gcm(modname,abort_message,42)
+        endif
+        firstcall=.false.
+      ENDIF
 
-
+      if (friction_type.eq.0) then
 c   calcul des composantes au carre du vent naturel
       do j=1,jjp1
@@ -96,5 +123,13 @@
          vcov(iip1,j,1)=vcov(1,j,1)
       enddo
+      endif ! of if (friction_type.eq.0)
 
+      if (friction_type.eq.1) then
+        do l=1,llm
+          ucov(:,:,l)=ucov(:,:,l)*(1.-pdt*kfrict(l))
+          vcov(:,:,l)=vcov(:,:,l)*(1.-pdt*kfrict(l))
+        enddo
+      endif
+      
       RETURN
       END
Index: /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/iniacademic.F
===================================================================
--- /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/iniacademic.F	(revision 1436)
+++ /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/iniacademic.F	(revision 1437)
@@ -9,4 +9,11 @@
       USE infotrac, ONLY : nqtot
       USE control_mod
+#ifdef CPP_IOIPSL
+      USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      USE ioipsl_getincom
+#endif
+      USE Write_Field
  
 
@@ -70,8 +77,14 @@
       REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
       REAL phi(ip1jmp1,llm)                  ! geopotentiel
-      REAL ddsin,tetarappelj,tetarappell,zsig
+      REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
       real tetajl(jjp1,llm)
       INTEGER i,j,l,lsup,ij
 
+      REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
+      REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
+      LOGICAL ok_geost             ! Initialisation vent geost. ou nul
+      LOGICAL ok_pv                ! Polar Vortex
+      REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex 
+      
       real zz,ran1
       integer idum
@@ -84,21 +97,16 @@
       if (planet_type=="earth") then
 c
+        ! initialize planet radius, rotation rate,...
+        call conf_planete
+
         time_0=0.
-        day_ref=0
+        day_ref=1
         annee_ref=0
 
         im         = iim
         jm         = jjm
-        day_ini    = 0
-        omeg       = 4.*asin(1.)/86400.
-        rad    = 6371229.
-        g      = 9.8
-        daysec = 86400.
+        day_ini    = 1
         dtvr    = daysec/REAL(day_step)
         zdtvr=dtvr
-        kappa  = 0.2857143
-        cpp    = 1004.70885
-        preff     = 101325.
-        pa        =  50000.
         etot0      = 0.
         ptot0      = 0.
@@ -129,30 +137,76 @@
         if (iflag_phys.eq.2) then
           ! initializations for the academic case
-          ps(:)=1.e5
-          phis(:)=0.
-c---------------------------------------------------------------------
-
-          taurappel=10.*daysec
-
-c---------------------------------------------------------------------
-c   Calcul de la temperature potentielle :
-c   --------------------------------------
-
+          
+          ! 1. local parameters
+          ! by convention, winter is in the southern hemisphere
+          ! Geostrophic wind or no wind?
+          ok_geost=.TRUE.
+          CALL getin('ok_geost',ok_geost)
+          ! Constants for Newtonian relaxation and friction
+          k_f=1.                !friction 
+          CALL getin('k_j',k_f)
+          k_f=1./(daysec*k_f)
+          k_c_s=4.  !cooling surface
+          CALL getin('k_c_s',k_c_s)
+          k_c_s=1./(daysec*k_c_s)
+          k_c_a=40. !cooling free atm
+          CALL getin('k_c_a',k_c_a)
+          k_c_a=1./(daysec*k_c_a)
+          ! Constants for Teta equilibrium profile
+          teta0=315.     ! mean Teta (S.H. 315K)
+          CALL getin('teta0',teta0)
+          ttp=200.       ! Tropopause temperature (S.H. 200K)
+          CALL getin('ttp',ttp)
+          eps=0.         ! Deviation to N-S symmetry(~0-20K)
+          CALL getin('eps',eps)
+          delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
+          CALL getin('delt_y',delt_y)
+          delt_z=10.     ! Vertical Gradient (S.H. 10K)
+          CALL getin('delt_z',delt_z)
+          ! Polar vortex
+          ok_pv=.false.
+          CALL getin('ok_pv',ok_pv)
+          phi_pv=-50.            ! Latitude of edge of vortex
+          CALL getin('phi_pv',phi_pv)
+          phi_pv=phi_pv*pi/180.
+          dphi_pv=5.             ! Width of the edge
+          CALL getin('dphi_pv',dphi_pv)
+          dphi_pv=dphi_pv*pi/180.
+          gam_pv=4.              ! -dT/dz vortex (in K/km)
+          CALL getin('gam_pv',gam_pv)
+          
+          ! 2. Initialize fields towards which to relax
+          ! Friction
+          knewt_g=k_c_a
           DO l=1,llm
-            zsig=ap(l)/preff+bp(l)
-            if (zsig.gt.0.3) then
-             lsup=l
-             tetarappell=1./8.*(-log(zsig)-.5)
-             DO j=1,jjp1
-             ddsin=sin(rlatu(j))-sin(pi/20.)
-             tetajl(j,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+tetarappell)
-             ENDDO
-            else
-c   Choix isotherme au-dessus de 300 mbar
-             do j=1,jjp1
-               tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa
-             enddo
-            endif ! of if (zsig.gt.0.3)
+           zsig=presnivs(l)/preff
+           knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
+           kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
+          ENDDO
+          DO j=1,jjp1
+            clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
+          ENDDO
+          
+          ! Potential temperature 
+          DO l=1,llm
+           zsig=presnivs(l)/preff
+           tetastrat=ttp*zsig**(-kappa)
+           tetapv=tetastrat
+           IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
+               tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
+           ENDIF
+           DO j=1,jjp1
+             ! Troposphere
+             ddsin=sin(rlatu(j))
+             tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
+     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
+             ! Profil stratospherique isotherme (+vortex)
+             w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
+             tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
+             tetajl(j,l)=MAX(tetajl(j,l),tetastrat)  
+           ENDDO
           ENDDO ! of DO l=1,llm
+ 
+!          CALL writefield('theta_eq',tetajl)
 
           do l=1,llm
@@ -165,37 +219,51 @@
           enddo
 
-c       call dump2d(jjp1,llm,tetajl,'TEQ   ')
-
-          CALL pression ( ip1jmp1, ap, bp, ps, p       )
-          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-          CALL massdair(p,masse)
-
-c  intialisation du vent et de la temperature
-          teta(:,:)=tetarappel(:,:)
-          CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
-          call ugeostr(phi,ucov)
-          vcov=0.
-          q(:,:,1   )=1.e-10
-          q(:,:,2   )=1.e-15
-          q(:,:,3:nqtot)=0.
-
-
-c   perturbation aleatoire sur la temperature
-          idum  = -1
-          zz = ran1(idum)
-          idum  = 0
-          do l=1,llm
-            do ij=iip2,ip1jm
-              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
-            enddo
-          enddo
-
-          do l=1,llm
-            do ij=1,ip1jmp1,iip1
-              teta(ij+iim,l)=teta(ij,l)
-            enddo
-          enddo
-
-
+          ! 3. Initialize fields (if necessary)
+          IF (.NOT. read_start) THEN
+            ! surface pressure
+            ps(:)=preff
+            ! ground geopotential
+            phis(:)=0.
+            
+            CALL pression ( ip1jmp1, ap, bp, ps, p       )
+            CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+            CALL massdair(p,masse)
+
+            ! bulk initialization of temperature
+            teta(:,:)=tetarappel(:,:)
+            
+            ! geopotential
+            CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
+            
+            ! winds
+            if (ok_geost) then
+              call ugeostr(phi,ucov)
+            else
+              ucov(:,:)=0.
+            endif
+            vcov(:,:)=0.
+            
+            ! bulk initialization of tracers
+            do i=1,nqtot
+              if (i.eq.1) q(:,:,i)=1.e-10
+              if (i.eq.2) q(:,:,i)=1.e-15
+              if (i.gt.2) q(:,:,i)=0.
+            enddo
+
+            ! add random perturbation to temperature
+            idum  = -1
+            zz = ran1(idum)
+            idum  = 0
+            do l=1,llm
+              do ij=iip2,ip1jm
+                teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
+              enddo
+            enddo
+
+            do l=1,llm
+              do ij=1,ip1jmp1,iip1
+                teta(ij+iim,l)=teta(ij,l)
+              enddo
+            enddo
 
 c     PRINT *,' Appel test_period avec tetarappel '
@@ -204,9 +272,11 @@
 c     CALL  test_period ( ucov,vcov,teta,q,p,phis )
 
-c   initialisation d'un traceur sur une colonne
-          j=jjp1*3/4
-          i=iip1/2
-          ij=(j-1)*iip1+i
-          q(ij,:,3)=1.
+           ! initialize a traceur on one column
+!          j=jjp1*3/4
+!          i=iip1/2
+!          ij=(j-1)*iip1+i
+!          q(ij,:,3)=1.
+
+          ENDIF ! of IF (.NOT. read_start)
         endif ! of if (iflag_phys.eq.2)
         
Index: /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/leapfrog.F
===================================================================
--- /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/leapfrog.F	(revision 1436)
+++ /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/leapfrog.F	(revision 1437)
@@ -430,10 +430,27 @@
 
       IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
-c   Calcul academique de la physique = Rappel Newtonien + friction 
-c   --------------------------------------------------------------
-       teta(:,:)=teta(:,:)
-     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
-       call friction(ucov,vcov,iphysiq*dtvr)
-      ENDIF
+!   Academic case : Simple friction and Newtonan relaxation 
+!   -------------------------------------------------------
+        DO l=1,llm   
+          DO ij=1,ip1jmp1
+           teta(ij,l)=teta(ij,l)-dtvr*
+     &      (teta(ij,l)-tetarappel(ij,l))*(knewt_g+knewt_t(l)*clat4(ij))
+          ENDDO
+        ENDDO ! of DO l=1,llm 
+          call friction(ucov,vcov,dtvr)
+        
+        ! Sponge layer (if any)
+        IF (ok_strato) THEN
+          dufi(:,:)=0.
+          dvfi(:,:)=0.
+          dtetafi(:,:)=0.
+          dqfi(:,:,:)=0.
+          dpfi(:)=0.
+          CALL top_bound(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
+          CALL addfi( dtvr, leapf, forward   ,
+     $                  ucov, vcov, teta , q   ,ps ,
+     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
+        ENDIF ! of IF (ok_strato) 
+      ENDIF ! of IF (iflag_phys.EQ.2)
 
 
Index: /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/academic.h
===================================================================
--- /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/academic.h	(revision 1436)
+++ /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/academic.h	(revision 1437)
@@ -1,5 +1,9 @@
 !
-! $Header$
+! $Id$
 !
-      real tetarappel(ip1jmp1,llm),taurappel
-      common/academic/tetarappel,taurappel
+      common/academic/tetarappel,knewt_t,kfrict,knewt_g,clat4
+      real :: tetarappel(ip1jmp1,llm)
+      real :: knewt_t(llm)
+      real :: kfrict(llm)
+      real :: knewt_g
+      real :: clat4(ip1jmp1)
Index: /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/comconst.h
===================================================================
--- /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/comconst.h	(revision 1436)
+++ /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/comconst.h	(revision 1437)
@@ -5,18 +5,33 @@
 ! INCLUDE comconst.h
 
-      COMMON/comconst/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,           &
-     & dtvr,daysec,                                                     &
+      COMMON/comconsti/im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl,          &
+     &                 iflag_top_bound
+      COMMON/comconstr/dtvr,daysec,                                     &
      & pi,dtphys,dtdiss,rad,r,cpp,kappa,cotot,unsim,g,omeg              &
      &                   ,dissip_factz,dissip_deltaz,dissip_zref        &
-     &                   ,iflag_top_bound,tau_top_bound
+     &                   ,tau_top_bound,                                &
+     & daylen,year_day,molmass
 
 
       INTEGER im,jm,lllm,imp1,jmp1,lllmm1,lllmp1,lcl
-      REAL dtvr,daysec
-      REAL pi,dtphys,dtdiss,rad,r,cpp,kappa
-      REAL cotot,unsim,g,omeg
+      REAL dtvr ! dynamical time step (in s)
+      REAL daysec !length (in s) of a standard day
+      REAL pi    ! something like 3.14159....
+      REAL dtphys ! (s) time step for the physics
+      REAL dtdiss ! (s) time step for the dissipation
+      REAL rad ! (m) radius of the planet
+      REAL r ! Gas constant R=8.31 J.K-1.mol-1
+      REAL cpp   ! Cp
+      REAL kappa ! kappa=R/Cp 
+      REAL cotot
+      REAL unsim ! = 1./iim
+      REAL g ! (m/s2) gravity
+      REAL omeg ! (rad/s) rotation rate of the planet
       REAL dissip_factz,dissip_deltaz,dissip_zref
       INTEGER iflag_top_bound
       REAL tau_top_bound
+      REAL daylen ! length of solar day, in 'standard' day length
+      REAL year_day ! Number of standard days in a year
+      REAL molmass ! (g/mol) molar mass of the atmosphere
 
 
Index: /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/conf_planete.F90
===================================================================
--- /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/conf_planete.F90	(revision 1437)
+++ /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/conf_planete.F90	(revision 1437)
@@ -0,0 +1,66 @@
+!
+! $Id$
+!
+SUBROUTINE conf_planete
+!
+#ifdef CPP_IOIPSL
+USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+USE ioipsl_getincom
+#endif
+IMPLICIT NONE
+!
+!
+!   Declarations :
+!   --------------
+#include "dimensions.h"
+#include "comconst.h"
+#include "comvert.h"
+!
+!   local:
+!   ------
+
+! ---------------------------------------------
+! Initialisations de constantes de la dynamique
+! ---------------------------------------------
+! Pi
+pi=2.*asin(1.)
+
+!Reference surface pressure (Pa)
+preff=101325.
+CALL getin('preff', preff)
+! Reference pressure at which hybrid coord. become purely pressure
+! pa=50000.
+pa=preff/2.
+CALL getin('pa', pa)
+! Gravity
+g=9.80665
+CALL getin('g',g)
+! Molar mass of the atmosphere
+molmass = 28.9644
+CALL getin('molmass',molmass)
+! kappa=R/Cp et Cp      
+kappa = 2./7.
+CALL getin('kappa',kappa)
+cpp=8.3145/molmass/kappa*1000.
+CALL getin('cpp',cpp)
+! Radius of the planet
+rad = 6371229. 
+CALL getin('radius',rad)
+! Length of a standard day (s)
+daysec=86400.
+CALL getin('daysec',daysec)
+! Rotation rate of the planet:
+! Length of a solar day, in standard days
+daylen = 1.
+CALL getin('daylen',daylen)
+! Number of days (standard) per year:
+year_day = 365.25
+CALL getin('year_day',year_day)
+! Omega
+! omeg=2.*pi/86400.
+omeg=2.*pi/daysec*(1./daylen+1./year_day)
+CALL getin('omeg',omeg)
+
+END SUBROUTINE conf_planete
Index: /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/friction_p.F
===================================================================
--- /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/friction_p.F	(revision 1436)
+++ /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/friction_p.F	(revision 1437)
@@ -6,17 +6,23 @@
       USE parallel
       USE control_mod
+#ifdef CPP_IOIPSL
+      USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      USE ioipsl_getincom
+#endif
       IMPLICIT NONE
 
-c=======================================================================
-c
-c
-c   Objet:
-c   ------
-c
-c  ***********
-c    Friction
-c  ***********
-c
-c=======================================================================
+!=======================================================================
+!
+!   Friction for the Newtonian case:
+!   --------------------------------
+!    2 possibilities (depending on flag 'friction_type'
+!     friction_type=0 : A friction that is only applied to the lowermost
+!                       atmospheric layer
+!     friction_type=1 : Friction applied on all atmospheric layer (but
+!       (default)       with stronger magnitude near the surface; see
+!                       iniacademic.F)
+!=======================================================================
 
 #include "dimensions.h"
@@ -24,16 +30,40 @@
 #include "comgeom2.h"
 #include "comconst.h"
-
-      REAL pdt
+#include "iniprint.h"
+#include "academic.h"
+
+! arguments:
+      REAL,INTENT(out) :: ucov( iip1,jjp1,llm )
+      REAL,INTENT(out) :: vcov( iip1,jjm,llm )
+      REAL,INTENT(in) :: pdt ! time step
+
+! local variables:
       REAL modv(iip1,jjp1),zco,zsi
       REAL vpn,vps,upoln,upols,vpols,vpoln
       REAL u2(iip1,jjp1),v2(iip1,jjm)
-      REAL ucov( iip1,jjp1,llm ),vcov( iip1,jjm,llm )
-      INTEGER  i,j
-      REAL cfric
-      parameter (cfric=1.e-5)
+      INTEGER  i,j,l
+      REAL,PARAMETER :: cfric=1.e-5
+      LOGICAL,SAVE :: firstcall=.true.
+      INTEGER,SAVE :: friction_type=1
+      CHARACTER(len=20) :: modname="friction_p"
+      CHARACTER(len=80) :: abort_message
+!$OMP THREADPRIVATE(firstcall,friction_type)
       integer :: jjb,jje
 
-
+!$OMP SINGLE
+      IF (firstcall) THEN
+        ! set friction type
+        call getin("friction_type",friction_type)
+        if ((friction_type.lt.0).or.(friction_type.gt.1)) then
+          abort_message="wrong friction type"
+          write(lunout,*)'Friction: wrong friction type',friction_type
+          call abort_gcm(modname,abort_message,42)
+        endif
+        firstcall=.false.
+      ENDIF
+!$OMP END SINGLE COPYPRIVATE(friction_type,firstcall)
+
+      if (friction_type.eq.0) then ! friction on first layer only
+!$OMP SINGLE
 c   calcul des composantes au carre du vent naturel
       jjb=jj_begin
@@ -138,4 +168,33 @@
          vcov(iip1,j,1)=vcov(1,j,1)
       enddo
+!$OMP END SINGLE
+      endif ! of if (friction_type.eq.0)
+
+      if (friction_type.eq.1) then
+       ! for ucov() 
+        jjb=jj_begin
+        jje=jj_end
+        if (pole_nord) jjb=jj_begin+1
+        if (pole_sud) jje=jj_end-1
+
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        do l=1,llm
+          ucov(1:iip1,jjb:jje,l)=ucov(1:iip1,jjb:jje,l)*
+     &                                  (1.-pdt*kfrict(l))
+        enddo
+!$OMP END DO NOWAIT
+        
+       ! for vcoc()
+        jjb=jj_begin
+        jje=jj_end
+        if (pole_sud) jje=jj_end-1
+        
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+        do l=1,llm
+          vcov(1:iip1,jjb:jje,l)=vcov(1:iip1,jjb:jje,l)*
+     &                                  (1.-pdt*kfrict(l))
+        enddo
+!$OMP END DO
+      endif ! of if (friction_type.eq.1)
 
       RETURN
Index: /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/iniacademic.F
===================================================================
--- /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/iniacademic.F	(revision 1436)
+++ /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/iniacademic.F	(revision 1437)
@@ -9,4 +9,11 @@
       USE infotrac, ONLY : nqtot
       USE control_mod
+#ifdef CPP_IOIPSL
+      USE IOIPSL
+#else
+! if not using IOIPSL, we still need to use (a local version of) getin
+      USE ioipsl_getincom
+#endif
+      USE Write_Field
  
 
@@ -70,8 +77,14 @@
       REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
       REAL phi(ip1jmp1,llm)                  ! geopotentiel
-      REAL ddsin,tetarappelj,tetarappell,zsig
+      REAL ddsin,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
       real tetajl(jjp1,llm)
       INTEGER i,j,l,lsup,ij
 
+      REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
+      REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
+      LOGICAL ok_geost             ! Initialisation vent geost. ou nul
+      LOGICAL ok_pv                ! Polar Vortex
+      REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex 
+      
       real zz,ran1
       integer idum
@@ -84,21 +97,16 @@
       if (planet_type=="earth") then
 c
+        ! initialize planet radius, rotation rate,...
+        call conf_planete
+
         time_0=0.
-        day_ref=0
+        day_ref=1
         annee_ref=0
 
         im         = iim
         jm         = jjm
-        day_ini    = 0
-        omeg       = 4.*asin(1.)/86400.
-        rad    = 6371229.
-        g      = 9.8
-        daysec = 86400.
+        day_ini    = 1
         dtvr    = daysec/REAL(day_step)
         zdtvr=dtvr
-        kappa  = 0.2857143
-        cpp    = 1004.70885
-        preff     = 101325.
-        pa        =  50000.
         etot0      = 0.
         ptot0      = 0.
@@ -129,30 +137,76 @@
         if (iflag_phys.eq.2) then
           ! initializations for the academic case
-          ps(:)=1.e5
-          phis(:)=0.
-c---------------------------------------------------------------------
-
-          taurappel=10.*daysec
-
-c---------------------------------------------------------------------
-c   Calcul de la temperature potentielle :
-c   --------------------------------------
-
+          
+          ! 1. local parameters
+          ! by convention, winter is in the southern hemisphere
+          ! Geostrophic wind or no wind?
+          ok_geost=.TRUE.
+          CALL getin('ok_geost',ok_geost)
+          ! Constants for Newtonian relaxation and friction
+          k_f=1.                !friction 
+          CALL getin('k_j',k_f)
+          k_f=1./(daysec*k_f)
+          k_c_s=4.  !cooling surface
+          CALL getin('k_c_s',k_c_s)
+          k_c_s=1./(daysec*k_c_s)
+          k_c_a=40. !cooling free atm
+          CALL getin('k_c_a',k_c_a)
+          k_c_a=1./(daysec*k_c_a)
+          ! Constants for Teta equilibrium profile
+          teta0=315.     ! mean Teta (S.H. 315K)
+          CALL getin('teta0',teta0)
+          ttp=200.       ! Tropopause temperature (S.H. 200K)
+          CALL getin('ttp',ttp)
+          eps=0.         ! Deviation to N-S symmetry(~0-20K)
+          CALL getin('eps',eps)
+          delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
+          CALL getin('delt_y',delt_y)
+          delt_z=10.     ! Vertical Gradient (S.H. 10K)
+          CALL getin('delt_z',delt_z)
+          ! Polar vortex
+          ok_pv=.false.
+          CALL getin('ok_pv',ok_pv)
+          phi_pv=-50.            ! Latitude of edge of vortex
+          CALL getin('phi_pv',phi_pv)
+          phi_pv=phi_pv*pi/180.
+          dphi_pv=5.             ! Width of the edge
+          CALL getin('dphi_pv',dphi_pv)
+          dphi_pv=dphi_pv*pi/180.
+          gam_pv=4.              ! -dT/dz vortex (in K/km)
+          CALL getin('gam_pv',gam_pv)
+          
+          ! 2. Initialize fields towards which to relax
+          ! Friction
+          knewt_g=k_c_a
           DO l=1,llm
-            zsig=ap(l)/preff+bp(l)
-            if (zsig.gt.0.3) then
-             lsup=l
-             tetarappell=1./8.*(-log(zsig)-.5)
-             DO j=1,jjp1
-             ddsin=sin(rlatu(j))-sin(pi/20.)
-             tetajl(j,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+tetarappell)
-             ENDDO
-            else
-c   Choix isotherme au-dessus de 300 mbar
-             do j=1,jjp1
-               tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa
-             enddo
-            endif ! of if (zsig.gt.0.3)
+           zsig=presnivs(l)/preff
+           knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
+           kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
+          ENDDO
+          DO j=1,jjp1
+            clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
+          ENDDO
+          
+          ! Potential temperature 
+          DO l=1,llm
+           zsig=presnivs(l)/preff
+           tetastrat=ttp*zsig**(-kappa)
+           tetapv=tetastrat
+           IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
+               tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
+           ENDIF
+           DO j=1,jjp1
+             ! Troposphere
+             ddsin=sin(rlatu(j))
+             tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
+     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
+             ! Profil stratospherique isotherme (+vortex)
+             w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
+             tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
+             tetajl(j,l)=MAX(tetajl(j,l),tetastrat)  
+           ENDDO
           ENDDO ! of DO l=1,llm
+ 
+!          CALL writefield('theta_eq',tetajl)
 
           do l=1,llm
@@ -165,37 +219,51 @@
           enddo
 
-c       call dump2d(jjp1,llm,tetajl,'TEQ   ')
-
-          CALL pression ( ip1jmp1, ap, bp, ps, p       )
-          CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-          CALL massdair(p,masse)
-
-c  intialisation du vent et de la temperature
-          teta(:,:)=tetarappel(:,:)
-          CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
-          call ugeostr(phi,ucov)
-          vcov=0.
-          q(:,:,1   )=1.e-10
-          q(:,:,2   )=1.e-15
-          q(:,:,3:nqtot)=0.
-
-
-c   perturbation aleatoire sur la temperature
-          idum  = -1
-          zz = ran1(idum)
-          idum  = 0
-          do l=1,llm
-            do ij=iip2,ip1jm
-              teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
-            enddo
-          enddo
-
-          do l=1,llm
-            do ij=1,ip1jmp1,iip1
-              teta(ij+iim,l)=teta(ij,l)
-            enddo
-          enddo
-
-
+          ! 3. Initialize fields (if necessary)
+          IF (.NOT. read_start) THEN
+            ! surface pressure
+            ps(:)=preff
+            ! ground geopotential
+            phis(:)=0.
+            
+            CALL pression ( ip1jmp1, ap, bp, ps, p       )
+            CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
+            CALL massdair(p,masse)
+
+            ! bulk initialization of temperature
+            teta(:,:)=tetarappel(:,:)
+            
+            ! geopotential
+            CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
+            
+            ! winds
+            if (ok_geost) then
+              call ugeostr(phi,ucov)
+            else
+              ucov(:,:)=0.
+            endif
+            vcov(:,:)=0.
+            
+            ! bulk initialization of tracers
+            do i=1,nqtot
+              if (i.eq.1) q(:,:,i)=1.e-10
+              if (i.eq.2) q(:,:,i)=1.e-15
+              if (i.gt.2) q(:,:,i)=0.
+            enddo
+
+            ! add random perturbation to temperature
+            idum  = -1
+            zz = ran1(idum)
+            idum  = 0
+            do l=1,llm
+              do ij=iip2,ip1jm
+                teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
+              enddo
+            enddo
+
+            do l=1,llm
+              do ij=1,ip1jmp1,iip1
+                teta(ij+iim,l)=teta(ij,l)
+              enddo
+            enddo
 
 c     PRINT *,' Appel test_period avec tetarappel '
@@ -204,9 +272,11 @@
 c     CALL  test_period ( ucov,vcov,teta,q,p,phis )
 
-c   initialisation d'un traceur sur une colonne
-          j=jjp1*3/4
-          i=iip1/2
-          ij=(j-1)*iip1+i
-          q(ij,:,3)=1.
+           ! initialize a traceur on one column
+!          j=jjp1*3/4
+!          i=iip1/2
+!          ij=(j-1)*iip1+i
+!          q(ij,:,3)=1.
+
+          ENDIF ! of IF (.NOT. read_start)
         endif ! of if (iflag_phys.eq.2)
         
Index: /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/leapfrog_p.F	(revision 1436)
+++ /LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/leapfrog_p.F	(revision 1437)
@@ -968,15 +968,14 @@
 
       IF(iflag_phys.EQ.2) THEN ! "Newtonian" case
-c   Calcul academique de la physique = Rappel Newtonien + fritcion 
-c   --------------------------------------------------------------
-cym       teta(:,:)=teta(:,:)
-cym     s  -iphysiq*dtvr*(teta(:,:)-tetarappel(:,:))/taurappel
+!   Academic case : Simple friction and Newtonan relaxation 
+!   -------------------------------------------------------
        ijb=ij_begin
        ije=ij_end
 !$OMP DO SCHEDULE(STATIC,OMP_CHUNK) 
        do l=1,llm
-       teta(ijb:ije,l)=teta(ijb:ije,l)
-     &  -iphysiq*dtvr*(teta(ijb:ije,l)-tetarappel(ijb:ije,l))/taurappel
-       enddo
+        teta(ijb:ije,l)=teta(ijb:ije,l)-dtvr*
+     &         (teta(ijb:ije,l)-tetarappel(ijb:ije,l))*
+     &                  (knewt_g+knewt_t(l)*clat4(ijb:ije))
+       enddo ! of do l=1,llm
 !$OMP END DO
 
@@ -987,8 +986,37 @@
        call WaitRequest(Request_Physic)     
 c$OMP BARRIER
-!$OMP MASTER
-       call friction_p(ucov,vcov,iphysiq*dtvr)
-!$OMP END MASTER
+       call friction_p(ucov,vcov,dtvr)
 !$OMP BARRIER
+
+        ! Sponge layer (if any)
+        IF (ok_strato) THEN
+          ! set dufi,dvfi,... to zero
+          ijb=ij_begin
+          ije=ij_end
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+          do l=1,llm
+            dufi(ijb:ije,l)=0
+            dtetafi(ijb:ije,l)=0
+            dqfi(ijb:ije,l,1:nqtot)=0
+          enddo
+!$OMP END DO
+!$OMP SINGLE
+          dpfi(ijb:ije)=0
+!$OMP END SINGLE
+          ijb=ij_begin
+          ije=ij_end
+          if (pole_sud) ije=ije-iip1
+!$OMP DO SCHEDULE(STATIC,OMP_CHUNK)
+          do l=1,llm
+            dvfi(ijb:ije,l)=0
+          enddo
+!$OMP END DO
+
+          CALL top_bound_p(vcov,ucov,teta,masse,dufi,dvfi,dtetafi)
+          CALL addfi_p( dtvr, leapf, forward   ,
+     $                  ucov, vcov, teta , q   ,ps ,
+     $                 dufi, dvfi, dtetafi , dqfi ,dpfi  )
+!$OMP BARRIER
+        ENDIF ! of IF (ok_strato) 
       ENDIF ! of IF(iflag_phys.EQ.2)
 
