MODULE etat0_mod
    CHARACTER(len=255),SAVE :: etat0_type

CONTAINS
  
  SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_q)
    USE icosa
    USE disvert_mod
    USE etat0_jablonowsky06_mod, ONLY : etat0_jablonowsky06=>etat0
    USE etat0_academic_mod, ONLY : etat0_academic=>etat0  
    USE etat0_dcmip1_mod, ONLY : etat0_dcmip1=>etat0
    USE etat0_dcmip2_mod, ONLY : etat0_dcmip2=>etat0
    USE etat0_dcmip3_mod, ONLY : etat0_dcmip3=>etat0  
    USE etat0_dcmip4_mod, ONLY : etat0_dcmip4=>etat0  
    USE etat0_dcmip5_mod, ONLY : etat0_dcmip5=>etat0  
    USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0  
    USE dynetat0_gcm_mod, ONLY : dynetat0_start=>etat0 
    USE dynetat0_hz_mod,  ONLY : dynetat0_hz=>etat0 

    IMPLICIT NONE
    TYPE(t_field),POINTER :: f_ps(:)
    TYPE(t_field),POINTER :: f_mass(:)
    TYPE(t_field),POINTER :: f_phis(:)
    TYPE(t_field),POINTER :: f_theta_rhodz(:)
    TYPE(t_field),POINTER :: f_u(:)
    TYPE(t_field),POINTER :: f_q(:)
    REAL(rstd),POINTER :: ps(:), mass(:,:)
    LOGICAL :: init_mass
    INTEGER :: ind,i,j,ij,l

    ! most etat0 routines set ps and not mass
    ! in that case and if caldyn_eta == eta_lag
    ! the initial distribution of mass is taken to be the same
    ! as what the mass coordinate would dictate
    ! however if etat0_XXX defines mass then the flag init_mass must be set to .FALSE.
    ! otherwise mass will be overwritten 
    init_mass = (caldyn_eta == eta_lag)

    etat0_type='jablonowsky06'
    CALL getin("etat0",etat0_type)
    
    SELECT CASE (TRIM(etat0_type))
    CASE ('jablonowsky06')
       CALL etat0_jablonowsky06(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
    CASE ('academic')
       CALL etat0_academic(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
    CASE ('heldsz')
	print*,"heldsz test case"
       CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
    CASE ('dcmip1')
       CALL etat0_dcmip1(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
       CALL etat0_dcmip2(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
    CASE ('dcmip3')
       CALL etat0_dcmip3(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
     CASE ('dcmip4')
       CALL etat0_dcmip4(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
     CASE ('dcmip5')
       CALL etat0_dcmip5(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
     CASE ('readnf_start') 
	  print*,"readnf_start used"	
       CALL dynetat0_start(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
	CASE ('readnf_hz') 
	  print*,"readnf_hz used"
       CALL dynetat0_hz(f_ps,f_phis,f_theta_rhodz,f_u,f_q) 
   CASE DEFAULT
       PRINT*, 'Bad selector for variable etat0 <',etat0_type, &
            '> options are <jablonowsky06>, <academic>, <dcmip[1-4]> '
       STOP
    END SELECT

    IF(init_mass) THEN ! initialize mass distribution using ps
       !$OMP BARRIER
       DO ind=1,ndomain
          CALL swap_dimensions(ind)
          CALL swap_geometry(ind)
          mass=f_mass(ind); ps=f_ps(ind)
          CALL compute_rhodz(.TRUE., ps, mass)
       END DO
    END IF

  END SUBROUTINE etat0
         
END MODULE etat0_mod