Index: trunk/NEMO/OPA_SRC/istate.F90
===================================================================
--- trunk/NEMO/OPA_SRC/istate.F90	(revision 479)
+++ trunk/NEMO/OPA_SRC/istate.F90	(revision 508)
@@ -4,4 +4,14 @@
    !! Ocean state   :  initial state setting
    !!=====================================================================
+   !! History :   4.0  !  89-12  (P. Andrich)  Original code
+   !!             5.0  !  91-11  (G. Madec)  rewritting
+   !!             6.0  !  96-01  (G. Madec)  terrain following coordinates
+   !!             8.0  !  01-09  (M. Levy, M. Ben Jelloul)  istate_eel
+   !!             8.0  !  01-09  (M. Levy, M. Ben Jelloul)  istate_uvg
+   !!             9.0  !  03-08  (G. Madec)  F90: Free form, modules
+   !!             9.0  !  03-09  (G. Madec, C. Talandier)  add EEL R5
+   !!             9.0  !  04-05  (A. Koch-Larrouy)  istate_gyre 
+   !!             9.0  !  06-07  (S. Masson)  distributed restart using iom
+   !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
@@ -13,5 +23,4 @@
    !!   istate_uvg    : initial velocity in geostropic balance
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE oce             ! ocean dynamics and active tracers 
    USE dom_oce         ! ocean space and time domain 
@@ -19,5 +28,4 @@
    USE ldftra_oce      ! ocean active tracers: lateral physics
    USE zdf_oce         ! ocean vertical physics
-   USE in_out_manager  ! I/O manager
    USE phycst          ! physical constants
    USE wzvmod          ! verctical velocity               (wzv     routine)
@@ -26,10 +34,11 @@
    USE restart         ! ocean restart                   (rst_read routine)
    USE solisl          ! ???
-
+   USE in_out_manager  ! I/O manager
+   USE iom
+   
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC istate_init   ! routine called by step.F90
+   PUBLIC   istate_init   ! routine called by step.F90
 
    !! * Substitutions
@@ -37,7 +46,7 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
    !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -48,46 +57,14 @@
       !!                   ***  ROUTINE istate_init  ***
       !! 
-      !! ** Purpose :   Initialization of the dynamics and tracers.
-      !!
-      !! ** Method  :
-      !!
-      !! History :
-      !!   4.0  !  91-03  ()  Original code
-      !!        !  91-11  (G. Madec)
-      !!   9.0  !  03-09  (G. Madec)  F90: Free form, modules, orthogonality
-      !!----------------------------------------------------------------------
-      USE iom
-      !! * Local declarations
-      !CT INTEGER ::   inum
-      !!----------------------------------------------------------------------
-
-
-      ! Initialization to zero
-      ! ----------------------
-
-      !     before fields       !       now fields        !      after fields       !
-      ;   ub   (:,:,:) = 0.e0   ;   un   (:,:,:) = 0.e0   ;   ua   (:,:,:) = 0.e0
-      ;   vb   (:,:,:) = 0.e0   ;   vn   (:,:,:) = 0.e0   ;   va   (:,:,:) = 0.e0
-      ;                         ;   wn   (:,:,:) = 0.e0   ;
-      ;   rotb (:,:,:) = 0.e0   ;   rotn (:,:,:) = 0.e0   ;
-      ;   hdivb(:,:,:) = 0.e0   ;   hdivn(:,:,:) = 0.e0   ;
-
-      ;   tb   (:,:,:) = 0.e0   ;   tn   (:,:,:) = 0.e0   ;   ta   (:,:,:) = 0.e0
-      ;   sb   (:,:,:) = 0.e0   ;   sn   (:,:,:) = 0.e0   ;   sa   (:,:,:) = 0.e0
+      !! ** Purpose :   Initialization of the dynamics and tracer fields.
+      !!----------------------------------------------------------------------
+
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'istate_ini : Initialization of the dynamics and tracers'
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
 
       rhd  (:,:,:) = 0.e0
       rhop (:,:,:) = 0.e0
       rn2  (:,:,:) = 0.e0 
-
-#if defined key_dynspg_rl
-      ! rigid-lid formulation
-      bsfb(:,:) = 0.e0      ! before barotropic stream-function
-      bsfn(:,:) = 0.e0      ! now    barotropic stream-function
-      bsfd(:,:) = 0.e0      ! barotropic stream-function trend
-#endif
-      ! free surface formulation
-      sshb(:,:) = 0.e0      ! before sea-surface height
-      sshn(:,:) = 0.e0      ! now    sea-surface height
-
 
       IF( ln_rstart ) THEN                    ! Restart from a file
@@ -100,4 +77,11 @@
          neuler = 0                              ! Set time-step indicator at nit000 (euler forward)
          adatrj = 0._wp
+         !                                       ! Initialization of ocean to zero
+         !     before fields       !       now fields          
+         ;   ub   (:,:,:) = 0.e0   ;   un   (:,:,:) = 0.e0   
+         ;   vb   (:,:,:) = 0.e0   ;   vn   (:,:,:) = 0.e0    
+         ;   rotb (:,:,:) = 0.e0   ;   rotn (:,:,:) = 0.e0  
+         ;   hdivb(:,:,:) = 0.e0   ;   hdivn(:,:,:) = 0.e0  
+         !
          IF( cp_cfg == 'eel' ) THEN
             CALL istate_eel                      ! EEL   configuration : start from pre-defined
@@ -107,11 +91,5 @@
             !                                    !                       and salinity fields 
          ELSE
-         !                                       ! Other configurations: Initial temperature and salinity fields
-
-         !CT CALL iom_open ('ssh', inum) 
-         !CT CALL iom_get( inum, jpdom_local, 'sshb', sshb )     ! free surface formulation (ssh)
-         !CT sshn(:,:) = sshb(:,:)
-         !CT CALL iom_close (inum)
-
+            !                                    ! Other configurations: Initial temperature and salinity fields
 #if defined key_dtatem
             CALL dta_tem( nit000 )                  ! read 3D temperature data
@@ -120,5 +98,5 @@
 #else
             IF(lwp) WRITE(numout,*)                 ! analytical temperature profile
-            IF(lwp) WRITE(numout,*)' Temperature initialization using an analytic profile'
+            IF(lwp) WRITE(numout,*)'             Temperature initialization using an analytic profile'
             CALL istate_tem
 #endif
@@ -130,5 +108,5 @@
             ! No salinity data
             IF(lwp)WRITE(numout,*)                  ! analytical salinity profile
-            IF(lwp)WRITE(numout,*)' Salinity initialisation using a constant value'
+            IF(lwp)WRITE(numout,*)'             Salinity initialisation using a constant value'
             CALL istate_sal
 #endif
@@ -139,5 +117,5 @@
       !                                       ! -----------------
       CALL wzv( nit000 )                         ! from horizontal divergence
-
+      !
    END SUBROUTINE istate_init
 
@@ -153,14 +131,8 @@
       !!
       !! References :  Philander ???
-      !!
-      !! History :
-      !!   4.0  !  89-12  (P. Andrich)  Original code
-      !!   6.0  !  96-01  (G. Madec)  terrain following coordinates
-      !!   9.0  !  02-09  (G. Madec)  F90: Free form
-      !!----------------------------------------------------------------------
-      !! * Local declarations
+      !!----------------------------------------------------------------------
       INTEGER :: ji, jj, jk
       !!----------------------------------------------------------------------
-
+      !
       IF(lwp) WRITE(numout,*)
       IF(lwp) WRITE(numout,*) 'istate_tem : initial temperature profile'
@@ -181,5 +153,5 @@
          &                 1     , jpi   , 5     , 1     , jpk   ,   &
          &                 1     , 1.    , numout                  )
-
+      !
    END SUBROUTINE istate_tem
 
@@ -194,10 +166,5 @@
       !!              
       !! ** Action  :   Initialize sn and sb
-      !!
-      !! History :
-      !!   4.0  !  89-12  (P. Andrich)  Original code
-      !!   8.5  !  02-09  (G. Madec)  F90: Free form
-      !!----------------------------------------------------------------------
-      !! * Local declarations
+      !!----------------------------------------------------------------------
       REAL(wp) ::   zsal = 35.50_wp
       !!----------------------------------------------------------------------
@@ -224,30 +191,19 @@
       !!              - set velocity field including horizontal divergence
       !!                and relative vorticity fields
-      !!
-      !! History :
-      !!   8.0  !  01-09  (M. Levy, M. Ben Jelloul)  read file for EEL 2 & 6
-      !!   9.0  !  03-09  (G. Madec, C. Talandier)  EEL 5
-      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
-      !!----------------------------------------------------------------------
-      !! * Modules used
+      !!----------------------------------------------------------------------
       USE eosbn2     ! eq. of state, Brunt Vaisala frequency (eos     routine)
       USE divcur     ! hor. divergence & rel. vorticity      (div_cur routine)
       USE iom
  
-      !! * Local declarations
       INTEGER  ::   inum              ! temporary logical unit
       INTEGER  ::   ji, jj, jk        ! dummy loop indices
       INTEGER  ::   ijloc
-      REAL(wp) ::   &
-         zh1, zh2, zslope, zcst,   &  ! temporary scalars
-         zfcor
-      REAL(wp) ::   &
-         zt1  = 12._wp,            &  ! surface temperature value (EEL R5)
-         zt2  =  2._wp,            &  ! bottom  temperature value (EEL R5)
-         zsal = 35.5_wp,           &  ! constant salinity (EEL R2, R5 and R6)
-         zueel = 0.1_wp               ! constant uniform zonal velocity (EEL R5)
+      REAL(wp) ::   zh1, zh2, zslope, zcst, zfcor   ! temporary scalars
+      REAL(wp) ::   zt1  = 12._wp,               &  ! surface temperature value (EEL R5)
+         &          zt2  =  2._wp,               &  ! bottom  temperature value (EEL R5)
+         &          zsal = 35.5_wp,              &  ! constant salinity (EEL R2, R5 and R6)
+         &          zueel = 0.1_wp                  ! constant uniform zonal velocity (EEL R5)
 # if ! defined key_dynspg_rl
-      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   &
-         zssh                         ! initial ssh over the global domain
+      REAL(wp), DIMENSION(jpiglo,jpjglo) ::   zssh  ! initial ssh over the global domain
 # endif
       !!----------------------------------------------------------------------
@@ -389,17 +345,8 @@
       !! ** Method  : - set temprature field
       !!              - set salinity field
-      !!
-      !! ** History :                                     
-      !!      9.0  !  04-05  (A. Koch-Larrouy)  Original code 
-      !!----------------------------------------------------------------------
-      !! * Modules used
-      USE iom
-
-      !! * Local variables
-      INTEGER  ::   inum              ! temporary logical unit
-      INTEGER, PARAMETER ::   &
-         ntsinit = 0         ! (0/1) (analytical/input data files) T&S initialization
-
+      !!----------------------------------------------------------------------
       INTEGER :: ji, jj, jk  ! dummy loop indices
+      INTEGER            ::   inum          ! temporary logical unit
+      INTEGER, PARAMETER ::   ntsinit = 0   ! (0/1) (analytical/input data files) T&S initialization
       !!----------------------------------------------------------------------
 
@@ -469,7 +416,5 @@
       ENDIF
 
-     
    END SUBROUTINE istate_gyre
-
 
 
@@ -484,11 +429,5 @@
       !!      surface to the bottom.
       !!                 p=integral [ rau*g dz ]
-      !!
-      !! History :
-      !!   8.1  !  01-09  (M. Levy, M. Ben Jelloul)  Original code
-      !!   8.5  !  02-09  (G. Madec)  F90: Free form
-      !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
-      !!----------------------------------------------------------------------
-      !! * Modules used
+      !!----------------------------------------------------------------------
       USE eosbn2          ! eq. of state, Brunt Vaisala frequency (eos     routine)
       USE dynspg          ! surface pressure gradient             (dyn_spg routine)
@@ -496,12 +435,8 @@
       USE lbclnk          ! ocean lateral boundary condition (or mpp link)
 
-      !! * Local declarations
       INTEGER ::   ji, jj, jk        ! dummy loop indices
       INTEGER ::   indic             ! ???
-      REAL(wp) ::   &
-         zmsv, zphv, zmsu, zphu,  &  ! temporary scalars
-         zalfg
-      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   &
-         zprn                        ! workspace
+      REAL(wp) ::   zmsv, zphv, zmsu, zphu, zalfg     ! temporary scalars
+      REAL(wp), DIMENSION (jpi,jpj,jpk) ::   zprn     ! workspace
       !!----------------------------------------------------------------------
 
@@ -514,9 +449,8 @@
 
       zalfg = 0.5 * grav * rau0
-      ! Surface value
-      zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) )
-
-      ! Vertical integration from the surface
-      DO jk = 2, jpkm1
+      
+      zprn(:,:,1) = zalfg * fse3w(:,:,1) * ( 1 + rhd(:,:,1) )       ! Surface value
+
+      DO jk = 2, jpkm1                                              ! Vertical integration from the surface
          zprn(:,:,jk) = zprn(:,:,jk-1)   &
             &         + zalfg * fse3w(:,:,jk) * ( 2. + rhd(:,:,jk) + rhd(:,:,jk-1) )
@@ -525,5 +459,4 @@
       ! Compute geostrophic balance
       ! ---------------------------
-
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
@@ -571,15 +504,9 @@
       ! after initializing u and v, we need to calculate the initial streamfunction bsf.
       ! Otherwise, only the trend will be computed and the model will blow up (inconsistency).
-      
       ! to do that, we call dyn_spg with a special trick:
-      ! we fill ua and va with the velocities divided by dt,
-      ! and the streamfunction will be brought to the right
-      ! value assuming the velocities have been set up in
-      ! one time step.
-      ! we then set bsfd to zero (first guess for next step
-      ! is d(psi)/dt = 0.)
-
-      !  sets up s false trend to calculate the barotropic
-      !  streamfunction.
+      ! we fill ua and va with the velocities divided by dt, and the streamfunction will be brought to the
+      ! right value assuming the velocities have been set up in one time step.
+      ! we then set bsfd to zero (first guess for next step is d(psi)/dt = 0.)
+      !  sets up s false trend to calculate the barotropic streamfunction.
 
       ua(:,:,:) = ub(:,:,:) / rdt
@@ -612,5 +539,5 @@
       hdivb(:,:,:) = hdivn(:,:,:)       ! set the before to the now value
       rotb (:,:,:) = rotn (:,:,:)       ! set the before to the now value
-
+      !
    END SUBROUTINE istate_uvg
 
