Index: LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F90
===================================================================
--- LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F	(revision 1455)
+++ LMDZ5/branches/LMDZ5V2.0-dev/libf/dyn3d/iniacademic.F90	(revision 1472)
@@ -2,284 +2,251 @@
 ! $Id$
 !
-c
-c
-      SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
-
-      USE filtreg_mod
-      USE infotrac, ONLY : nqtot
-      USE control_mod, ONLY: day_step,planet_type
+SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
+
+  USE filtreg_mod
+  USE infotrac, ONLY : nqtot
+  USE control_mod, ONLY: day_step,planet_type
 #ifdef CPP_IOIPSL
-      USE IOIPSL
+  USE IOIPSL
 #else
-! if not using IOIPSL, we still need to use (a local version of) getin
-      USE ioipsl_getincom
+  ! if not using IOIPSL, we still need to use (a local version of) getin
+  USE ioipsl_getincom
 #endif
-      USE Write_Field
- 
-
-c%W%    %G%
-c=======================================================================
-c
-c   Author:    Frederic Hourdin      original: 15/01/93
-c   -------
-c
-c   Subject:
-c   ------
-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 "academic.h"
-#include "ener.h"
-#include "temps.h"
-#include "iniprint.h"
-#include "logic.h"
-
-c   Arguments:
-c   ----------
-
-      real time_0
-
-c   variables dynamiques
-      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
-      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
-      REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
-      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 phi(ip1jmp1,llm)                  ! geopotentiel
-      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
-
-      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
-
-c-----------------------------------------------------------------------
-! 1. Initializations for Earth-like case
-! --------------------------------------
-c
-        ! initialize planet radius, rotation rate,...
-        call conf_planete
-
-        time_0=0.
-        day_ref=1
-        annee_ref=0
-
-        im         = iim
-        jm         = jjm
-        day_ini    = 1
-        dtvr    = daysec/REAL(day_step)
-        zdtvr=dtvr
-        etot0      = 0.
-        ptot0      = 0.
-        ztot0      = 0.
-        stot0      = 0.
-        ang0       = 0.
-
-        if (llm.eq.1) then
-          ! specific initializations for the shallow water case
-          kappa=1
+  USE Write_Field
+
+  !   Author:    Frederic Hourdin      original: 15/01/93
+
+  IMPLICIT NONE
+
+  !   Declararations:
+  !   ---------------
+
+  include "dimensions.h"
+  include "paramet.h"
+  include "comvert.h"
+  include "comconst.h"
+  include "comgeom.h"
+  include "academic.h"
+  include "ener.h"
+  include "temps.h"
+  include "iniprint.h"
+  include "logic.h"
+
+  !   Arguments:
+  !   ----------
+
+  real time_0
+
+  !   variables dynamiques
+  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
+  REAL teta(ip1jmp1,llm)                 ! temperature potentielle
+  REAL q(ip1jmp1,llm,nqtot)               ! champs advectes
+  REAL ps(ip1jmp1)                       ! pression  au sol
+  REAL masse(ip1jmp1,llm)                ! masse d'air
+  REAL phis(ip1jmp1)                     ! geopotentiel au sol
+
+  !   Local:
+  !   ------
+
+  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 phi(ip1jmp1,llm)                  ! geopotentiel
+  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
+
+  REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr
+
+  !-----------------------------------------------------------------------
+  ! 1. Initializations for Earth-like case
+  ! --------------------------------------
+  !
+  ! initialize planet radius, rotation rate,...
+  call conf_planete
+
+  time_0=0.
+  day_ref=1
+  annee_ref=0
+
+  im         = iim
+  jm         = jjm
+  day_ini    = 1
+  dtvr    = daysec/REAL(day_step)
+  zdtvr=dtvr
+  etot0      = 0.
+  ptot0      = 0.
+  ztot0      = 0.
+  stot0      = 0.
+  ang0       = 0.
+
+  if (llm == 1) then
+     ! specific initializations for the shallow water case
+     kappa=1
+  endif
+
+  CALL iniconst
+  CALL inigeom
+  CALL inifilr
+
+  if (llm == 1) then
+     ! initialize fields for the shallow water case, if required
+     if (.not.read_start) then
+        phis(:)=0.
+        q(:,:,:)=0
+        CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
+     endif
+  endif
+
+  academic_case: if (iflag_phys == 2) then
+     ! initializations
+
+     ! 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=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
+
+     !          CALL writefield('theta_eq',tetajl)
+
+     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
+
+     ! 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
-        
-        CALL iniconst
-        CALL inigeom
-        CALL inifilr
-
-        if (llm.eq.1) then
-          ! initialize fields for the shallow water case, if required
-          if (.not.read_start) then
-            phis(:)=0.
-            q(:,:,:)=0
-            CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
-          endif
-        endif
-
-        if (iflag_phys.eq.2) then
-          ! initializations for the academic case
-          
-!         if (planet_type=="earth") then
-
-          ! 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=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
-            do j=1,jjp1
-              do i=1,iip1
-                 ij=(j-1)*iip1+i
-                 tetarappel(ij,l)=tetajl(j,l)
-              enddo
-            enddo
-          enddo
-
-
-!         else
-!          write(lunout,*)"iniacademic: planet types other than earth",
-!     &                   " not implemented (yet)."
-!          stop
-!         endif ! of if (planet_type=="earth")
-
-          ! 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
-            if (planet_type=="earth") then
-              ! Earth: first two tracers will be water
-              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
-            else
-              q(:,:,:)=0
-            endif ! of if (planet_type=="earth")
-
-            ! 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
-
-            ! maintain periodicity in longitude
-            do l=1,llm
-              do ij=1,ip1jmp1,iip1
-                teta(ij+iim,l)=teta(ij,l)
-              enddo
-            enddo
-
-          ENDIF ! of IF (.NOT. read_start)
-        endif ! of if (iflag_phys.eq.2)
-        
-      END
-c-----------------------------------------------------------------------
+        vcov(:,:)=0.
+
+        ! bulk initialization of tracers
+        if (planet_type=="earth") then
+           ! Earth: first two tracers will be water
+           do i=1,nqtot
+              if (i == 1) q(:,:,i)=1.e-10
+              if (i == 2) q(:,:,i)=1.e-15
+              if (i.gt.2) q(:,:,i)=0.
+           enddo
+        else
+           q(:,:,:)=0
+        endif ! of if (planet_type=="earth")
+
+        ! 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
+
+        ! maintain periodicity in longitude
+        do l=1,llm
+           do ij=1,ip1jmp1,iip1
+              teta(ij+iim,l)=teta(ij,l)
+           enddo
+        enddo
+
+     ENDIF ! of IF (.NOT. read_start)
+  endif academic_case
+
+END SUBROUTINE iniacademic
