Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/calfis_p.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/calfis_p.F	(revision 1357)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/calfis_p.F	(revision 1363)
@@ -108,4 +108,5 @@
 #include "comvert.h"
 #include "comgeom2.h"
+#include "iniprint.h"
 #ifdef CPP_MPI
       include 'mpif.h'
@@ -258,8 +259,9 @@
         debut = .TRUE.
         IF (ngridmx.NE.2+(jjm-1)*iim) THEN
-         PRINT*,'STOP dans calfis'
-         PRINT*,'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
-         PRINT*,'  ngridmx  jjm   iim   '
-         PRINT*,ngridmx,jjm,iim
+         write(lunout,*) 'STOP dans calfis'
+         write(lunout,*) 
+     &   'La dimension ngridmx doit etre egale a 2 + (jjm-1)*iim'
+         write(lunout,*) '  ngridmx  jjm   iim   '
+         write(lunout,*) ngridmx,jjm,iim
          STOP
         ENDIF
@@ -628,5 +630,7 @@
 #ifdef CPP_EARTH
 
-      print*,'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
+!$OMP MASTER
+      write(lunout,*) 'PHYSIQUE AVEC NSPLIT_PHYS=',nsplit_phys
+!$OMP END MASTER
       zdt_split=dtphys/nsplit_phys
       zdufic_omp(:,:)=0.
@@ -1107,5 +1111,6 @@
 
 #else
-      write(*,*) "calfis_p: for now can only work with parallel physics"
+      write(lunout,*)
+     & "calfis_p: for now can only work with parallel physics"
       stop
 #endif
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/exner_hyb.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/exner_hyb.F	(revision 1357)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/exner_hyb.F	(revision 1363)
@@ -1,4 +1,4 @@
 !
-! $Header$
+! $Id $
 !
       SUBROUTINE  exner_hyb ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
@@ -51,5 +51,33 @@
       REAL SSUM
 c
+
+      if (llm.eq.1) then
+        ! Specific behaviour for Shallow Water (1 vertical layer) case
       
+        ! Sanity checks
+        if (kappa.ne.1) then
+          call abort_gcm("exner_hyb",
+     &    "kappa!=1 , but running in Shallow Water mode!!",42)
+        endif
+        if (cpp.ne.r) then
+        call abort_gcm("exner_hyb",
+     &    "cpp!=r , but running in Shallow Water mode!!",42)
+        endif
+        
+        ! Compute pks(:),pk(:),pkf(:)
+        
+        DO   ij  = 1, ngrid
+          pks(ij) = (cpp/preff) * ps(ij) 
+          pk(ij,1) = .5*pks(ij)
+        ENDDO
+        
+        CALL SCOPY   ( ngrid * llm, pk, 1, pkf, 1 )
+        CALL filtreg ( pkf, jmp1, llm, 2, 1, .TRUE., 1 ) 
+        
+        ! our work is done, exit routine
+        return
+      endif ! of if (llm.eq.1)
+
+     
       unpl2k    = 1.+ 2.* kappa
 c
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/exner_hyb_p.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/exner_hyb_p.F	(revision 1357)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/exner_hyb_p.F	(revision 1363)
@@ -1,2 +1,5 @@
+!
+! $Id $
+!
       SUBROUTINE  exner_hyb_p ( ngrid, ps, p,alpha,beta, pks, pk, pkf )
 c
@@ -51,5 +54,67 @@
       INTEGER ije,ijb,jje,jjb
 c
-c$OMP BARRIER           
+c$OMP BARRIER
+
+      if (llm.eq.1) then
+        ! Specific behaviour for Shallow Water (1 vertical layer) case
+      
+        ! Sanity checks
+        if (kappa.ne.1) then
+          call abort_gcm("exner_hyb",
+     &    "kappa!=1 , but running in Shallow Water mode!!",42)
+        endif
+        if (cpp.ne.r) then
+        call abort_gcm("exner_hyb",
+     &    "cpp!=r , but running in Shallow Water mode!!",42)
+        endif
+        
+        ! Compute pks(:),pk(:),pkf(:)
+        ijb=ij_begin
+        ije=ij_end
+!$OMP DO SCHEDULE(STATIC)
+        DO ij=ijb, ije
+          pks(ij)=(cpp/preff)*ps(ij)
+          pk(ij,1) = .5*pks(ij)
+          pkf(ij,1)=pk(ij,1)
+        ENDDO
+!$OMP ENDDO
+
+!$OMP MASTER
+      if (pole_nord) then
+        DO  ij   = 1, iim
+          ppn(ij) = aire(   ij   ) * pks(  ij     )
+        ENDDO
+        xpn      = SSUM(iim,ppn,1) /apoln
+  
+        DO ij   = 1, iip1
+          pks(   ij     )  =  xpn
+          pk(ij,1) = .5*pks(ij)
+          pkf(ij,1)=pk(ij,1)
+        ENDDO
+      endif
+      
+      if (pole_sud) then
+        DO  ij   = 1, iim
+          pps(ij) = aire(ij+ip1jm) * pks(ij+ip1jm )
+        ENDDO
+        xps      = SSUM(iim,pps,1) /apols 
+  
+        DO ij   = 1, iip1
+          pks( ij+ip1jm )  =  xps
+          pk(ij+ip1jm,1)=.5*pks(ij+ip1jm)
+          pkf(ij+ip1jm,1)=pk(ij+ip1jm,1)
+        ENDDO
+      endif
+!$OMP END MASTER
+
+        jjb=jj_begin
+        jje=jj_end
+        CALL filtreg_p ( pkf,jjb,jje, jmp1, llm, 2, 1, .TRUE., 1 )
+
+        ! our work is done, exit routine
+        return
+      endif ! of if (llm.eq.1)
+
+
       unpl2k    = 1.+ 2.* kappa
 c
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/iniacademic.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/iniacademic.F	(revision 1357)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/iniacademic.F	(revision 1363)
@@ -8,4 +8,6 @@
       USE filtreg_mod
       USE infotrac, ONLY : nqtot
+      USE control_mod
+ 
 
 c%W%    %G%
@@ -31,5 +33,4 @@
 c
 c=======================================================================
-      USE control_mod
       IMPLICIT NONE
 c-----------------------------------------------------------------------
@@ -46,4 +47,5 @@
 #include "temps.h"
 #include "iniprint.h"
+#include "logic.h"
 
 c   Arguments:
@@ -55,5 +57,5 @@
       REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
       REAL teta(ip1jmp1,llm)                 ! temperature potentielle
-      REAL q(ip1jmp1,llm,nqtot)              ! champs advectes
+      REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
       REAL ps(ip1jmp1)                       ! pression  au sol
       REAL masse(ip1jmp1,llm)                ! masse d'air
@@ -84,5 +86,5 @@
         time_0=0.
         day_ref=0
-	annee_ref=0
+        annee_ref=0
 
         im         = iim
@@ -105,13 +107,31 @@
         ang0       = 0.
 
+        if (llm.eq.1) then
+          ! specific initializations for the shallow water case
+          kappa=1
+        endif
+        
         CALL iniconst
         CALL inigeom
         CALL inifilr
 
-        ps=0.
-        phis=0.
+        if (llm.eq.1) then
+          ! initialize fields for the shallow water case, if required
+          if (.not.read_start) then
+            phis(:)=0.
+            q(:,:,1)=1.e-10
+            q(:,:,2)=1.e-15
+            q(:,:,3:nqtot)=0.
+            CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
+          endif
+        endif
+
+        if (iflag_phys.eq.2) then
+          ! initializations for the academic case
+          ps(:)=1.e5
+          phis(:)=0.
 c---------------------------------------------------------------------
 
-        taurappel=10.*daysec
+          taurappel=10.*daysec
 
 c---------------------------------------------------------------------
@@ -119,63 +139,61 @@
 c   --------------------------------------
 
-        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
+          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
+             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)
-        ENDDO ! of DO l=1,llm
-
-        do l=1,llm
-           do j=1,jjp1
+             do j=1,jjp1
+               tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa
+             enddo
+            endif ! of if (zsig.gt.0.3)
+          ENDDO ! of DO l=1,llm
+
+          do l=1,llm
+            do j=1,jjp1
               do i=1,iip1
                  ij=(j-1)*iip1+i
                  tetarappel(ij,l)=tetajl(j,l)
               enddo
-           enddo
-        enddo
+            enddo
+          enddo
 
 c       call dump2d(jjp1,llm,tetajl,'TEQ   ')
 
-        ps=1.e5
-        phis=0.
-        CALL pression ( ip1jmp1, ap, bp, ps, p       )
-        CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
-        CALL massdair(p,masse)
+          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.
+          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
+          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
+            enddo
+          enddo
+
+          do l=1,llm
+            do ij=1,ip1jmp1,iip1
               teta(ij+iim,l)=teta(ij,l)
-           enddo
-        enddo
+            enddo
+          enddo
 
 
@@ -187,9 +205,10 @@
 
 c   initialisation d'un traceur sur une colonne
-        j=jjp1*3/4
-        i=iip1/2
-        ij=(j-1)*iip1+i
-        q(ij,:,3)=1.
-      
+          j=jjp1*3/4
+          i=iip1/2
+          ij=(j-1)*iip1+i
+          q(ij,:,3)=1.
+        endif ! of if (iflag_phys.eq.2)
+        
       else
         write(lunout,*)"iniacademic: planet types other than earth",
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/leapfrog_p.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/leapfrog_p.F	(revision 1357)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/leapfrog_p.F	(revision 1363)
@@ -1354,4 +1354,5 @@
                 call Gather_Field(ucov,ip1jmp1,llm,0)
                 call Gather_Field(teta,ip1jmp1,llm,0)
+                call Gather_Field(pk,ip1jmp1,llm,0)
 		call Gather_Field(phi,ip1jmp1,llm,0)
                 do iq=1,nqtot
@@ -1559,4 +1560,5 @@
                 call Gather_Field(ucov,ip1jmp1,llm,0)
                 call Gather_Field(teta,ip1jmp1,llm,0)
+                call Gather_Field(pk,ip1jmp1,llm,0)
 		call Gather_Field(phi,ip1jmp1,llm,0)
                 do iq=1,nqtot
@@ -1624,5 +1626,5 @@
 		call Gather_Field(phis,ip1jmp1,1,0)
                 if (mpi_rank==0) then
-	         CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
+                 CALL writehist(itau,vcov,ucov,teta,phi,q,masse,ps,phis)
 		endif
 !                CALL writehist_p(histid, histvid, itau,vcov , 
Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/sw_case_williamson91_6.F
===================================================================
--- LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/sw_case_williamson91_6.F	(revision 1363)
+++ LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3dpar/sw_case_williamson91_6.F	(revision 1363)
@@ -0,0 +1,140 @@
+!
+! $Id $
+!
+      SUBROUTINE sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
+
+c=======================================================================
+c
+c   Author:    Thomas Dubos      original: 26/01/2010
+c   -------
+c
+c   Subject:
+c   ------
+c   Realise le cas-test 6 de Williamson et al. (1991) : onde de Rossby-Haurwitz
+c
+c   Method:
+c   --------
+c
+c   Interface:
+c   ----------
+c
+c      Input:
+c      ------
+c
+c      Output:
+c      -------
+c
+c=======================================================================
+      IMPLICIT NONE
+c-----------------------------------------------------------------------
+c   Declararations:
+c   ---------------
+
+#include "dimensions.h"
+#include "paramet.h"
+#include "comvert.h"
+#include "comconst.h"
+#include "comgeom.h"
+#include "iniprint.h"
+
+c   Arguments:
+c   ----------
+
+c   variables dynamiques
+      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
+      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
+      REAL ps(ip1jmp1)                       ! pression  au sol
+      REAL masse(ip1jmp1,llm)                ! masse d'air
+      REAL phis(ip1jmp1)                     ! geopotentiel au sol
+
+c   Local:
+c   ------
+
+      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
+      REAL pks(ip1jmp1)                      ! exner au  sol
+      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
+      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
+      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
+
+      REAL :: sinth,costh,costh2, Ath,Bth,Cth, lon,dps
+      INTEGER i,j,ij
+
+      REAL, PARAMETER    :: rho=1 ! masse volumique de l'air (arbitraire)
+      REAL, PARAMETER    :: K    = 7.848e-6  ! K = \omega
+      REAL, PARAMETER    :: gh0  = 9.80616 * 8e3 
+      INTEGER, PARAMETER :: R0=4, R1=R0+1, R2=R0+2         ! mode 4
+c NB : rad = 6371220 dans W91 (6371229 dans LMDZ)
+c      omeg = 7.292e-5 dans W91 (7.2722e-5 dans LMDZ)
+ 
+      IF(0==0) THEN
+c Williamson et al. (1991) : onde de Rossby-Haurwitz
+         teta = preff/rho/cpp
+c geopotentiel (pression de surface)
+         do j=1,jjp1
+            costh2 = cos(rlatu(j))**2
+            Ath = (R0+1)*(costh2**2) + (2*R0*R0-R0-2)*costh2 - 2*R0*R0
+            Ath = .25*(K**2)*(costh2**(R0-1))*Ath
+            Ath = .5*K*(2*omeg+K)*costh2 + Ath 
+            Bth = (R1*R1+1)-R1*R1*costh2
+            Bth = 2*(omeg+K)*K/(R1*R2) * (costh2**(R0/2))*Bth
+            Cth = R1*costh2 - R2
+            Cth = .25*K*K*(costh2**R0)*Cth
+            do i=1,iip1
+               ij=(j-1)*iip1+i
+               lon = rlonv(i)
+               dps = Ath + Bth*cos(R0*lon) + Cth*cos(2*R0*lon)
+               ps(ij) = rho*(gh0 + (rad**2)*dps)
+            enddo
+         enddo
+         write(lunout,*) 'W91 ps', MAXVAL(ps), MINVAL(ps)
+c vitesse zonale ucov
+         do j=1,jjp1
+            costh  = cos(rlatu(j))
+            costh2 = costh**2
+            Ath = rad*K*costh
+            Bth = R0*(1-costh2)-costh2
+            Bth = rad*K*Bth*(costh**(R0-1))
+            do i=1,iip1
+               ij=(j-1)*iip1+i
+               lon = rlonu(i)
+               ucov(ij,1) = (Ath + Bth*cos(R0*lon))
+            enddo
+         enddo
+         write(lunout,*) 'W91 u', MAXVAL(ucov(:,1)), MINVAL(ucov(:,1))
+         ucov(:,1)=ucov(:,1)*cu
+c vitesse meridienne vcov
+         do j=1,jjm
+            sinth  = sin(rlatv(j))
+            costh  = cos(rlatv(j))
+            Ath = -rad*K*R0*sinth*(costh**(R0-1))
+            do i=1,iip1
+               ij=(j-1)*iip1+i
+               lon = rlonv(i)
+               vcov(ij,1) = Ath*sin(R0*lon)
+            enddo
+         enddo
+         write(lunout,*) 'W91 v', MAXVAL(vcov(:,1)), MINVAL(vcov(:,1))
+         vcov(:,1)=vcov(:,1)*cv
+         
+c         ucov=0
+c         vcov=0
+      ELSE
+c test non-tournant, onde se propageant en latitude
+         do j=1,jjp1
+            do i=1,iip1
+               ij=(j-1)*iip1+i
+               ps(ij) = 1e5*(1 + .1*exp(-100*(1+sin(rlatu(j)))**2) )
+            enddo
+         enddo
+         
+c     rho = preff/(cpp*teta)
+         teta = .01*preff/cpp   ! rho = 100 ; phi = ps/rho = 1e3 ; c=30 m/s = 2600 km/j = 23 degres / j
+         ucov=0.
+         vcov=0.
+      END IF      
+      
+      CALL pression ( ip1jmp1, ap, bp, ps, p       )
+      CALL massdair(p,masse)
+
+      END
+c-----------------------------------------------------------------------
