MODULE etat0_mod
  USE icosa
  USE etat0_collocated_mod
  IMPLICIT NONE 	
  PRIVATE

  PUBLIC :: etat0, init_etat0

! Important notes for OpenMP
! When etat0 is called, vertical OpenMP parallelism is deactivated. 
! Therefore only the omp_level_master thread must work, i.e. :
!   !$OMP BARRIER
!    DO ind=1,ndomain
!      IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
!      ...
!    END DO
!   !$OMP BARRIER
! There MUST be NO OMP BARRIER inside the DO-LOOP or any routine it calls.

CONTAINS
  
  SUBROUTINE init_etat0
    USE etat0_database_mod, ONLY: init_etat0_database => init_etat0 
    USE etat0_start_file_mod, ONLY: init_etat0_start_file => init_etat0 
    USE etat0_heldsz_mod, ONLY: init_etat0_held_suarez => init_etat0 
    
    CALL getin("etat0",etat0_type)

    SELECT CASE (TRIM(etat0_type))
      CASE ('isothermal')
      CASE ('temperature_profile')
      CASE ('jablonowsky06')
      CASE ('dcmip5')
      CASE ('williamson91.6')
      CASE ('start_file')
        CALL init_etat0_start_file
      CASE ('database')
        CALL init_etat0_database
      CASE ('academic')
      CASE ('held_suarez')
         CALL init_etat0_held_suarez
      CASE ('venus')
      CASE ('dcmip1')
      CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
      CASE ('dcmip3')
      CASE ('dcmip4')
      CASE ('dcmip2016_baroclinic_wave')
      CASE ('dcmip2016_cyclone')
      CASE ('dcmip2016_supercell')
      CASE ('bubble')
      CASE DEFAULT
         PRINT*, 'Bad selector for variable etat0 <',TRIM(etat0_type),'>'// &
            ' options are  <isothermal>, <temperature_profile>, <jablonowsky06>, <dcmip5>, <williamson91.6>,'& 
                         //' <start_file>, <database>, <academic>, <held_suarez>, <venus>, <dcmip1>,'         &
                         //' <dcmip2_mountain,dcmip2_schaer_noshear,dcmip2_schaer_shear>, <dcmip3>, <dcmip4>,'&
                         //' <dcmip2016_baroclinic_wave>, <dcmip2016_cyclone>, <dcmip2016_supercell>', 'bubble'
         STOP
    END SELECT

  END SUBROUTINE init_etat0

  SUBROUTINE etat0(f_ps,f_mass,f_phis,f_theta_rhodz,f_u, f_geopot,f_w, f_q)
    USE omp_para, ONLY : is_omp_level_master
    USE mpipara, ONLY : is_mpi_root
    USE disvert_mod, ONLY : caldyn_eta, eta_lag
    USE compute_diagnostics_mod, ONLY : compute_rhodz
    ! Generic interface
    USE etat0_isothermal_mod, ONLY : getin_etat0_isothermal=>getin_etat0
    USE etat0_dcmip1_mod, ONLY : getin_etat0_dcmip1=>getin_etat0
    USE etat0_dcmip2_mod, ONLY : getin_etat0_dcmip2=>getin_etat0
    USE etat0_dcmip4_mod, ONLY : getin_etat0_dcmip4=>getin_etat0
    USE etat0_dcmip5_mod, ONLY : getin_etat0_dcmip5=>getin_etat0
    USE etat0_bubble_mod, ONLY : getin_etat0_bubble=>getin_etat0
    USE etat0_williamson_mod, ONLY : getin_etat0_williamson=>getin_etat0
    USE etat0_temperature_mod, ONLY: getin_etat0_temperature=>getin_etat0
    USE etat0_dcmip2016_baroclinic_wave_mod, ONLY : getin_etat0_dcmip2016_baroclinic_wave=>getin_etat0
    USE etat0_dcmip2016_cyclone_mod, ONLY : getin_etat0_dcmip2016_cyclone=>getin_etat0
    USE etat0_dcmip2016_supercell_mod, ONLY : getin_etat0_dcmip2016_supercell=>getin_etat0
    ! Ad hoc interfaces
    USE etat0_academic_mod, ONLY : etat0_academic=>etat0
    USE etat0_heldsz_mod, ONLY : etat0_heldsz=>etat0
    USE etat0_venus_mod,  ONLY : etat0_venus=>etat0
    USE etat0_database_mod, ONLY : etat0_database=>etat0
    USE etat0_start_file_mod, ONLY : etat0_start_file=>etat0  

    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_geopot(:)
    TYPE(t_field),POINTER :: f_w(:)
    TYPE(t_field),POINTER :: f_q(:)
    
    REAL(rstd),POINTER :: ps(:), mass(:,:)
    LOGICAL :: autoinit_mass, autoinit_geopot, collocated
    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 autoinit_mass must be set to .FALSE.
    ! otherwise mass will be overwritten 
    autoinit_mass = (caldyn_eta == eta_lag)

    etat0_type='jablonowsky06'
    CALL getin("etat0",etat0_type)
    
    !------------------- Generic interface ---------------------
    collocated=.TRUE.
    SELECT CASE (TRIM(etat0_type))
    CASE ('isothermal')
       CALL getin_etat0_isothermal
    CASE ('temperature_profile')
       CALL getin_etat0_temperature
    CASE ('jablonowsky06')
    CASE ('dcmip1')
        CALL getin_etat0_dcmip1
    CASE ('dcmip2_mountain','dcmip2_schaer_noshear','dcmip2_schaer_shear')
       CALL getin_etat0_dcmip2
    CASE ('dcmip3')
    CASE ('dcmip4')
        CALL getin_etat0_dcmip4
    CASE ('dcmip5')
        CALL getin_etat0_dcmip5
    CASE ('bubble')
        CALL getin_etat0_bubble
    CASE ('williamson91.6')
       autoinit_mass=.FALSE.
       CALL getin_etat0_williamson
    CASE ('dcmip2016_baroclinic_wave')
        CALL getin_etat0_dcmip2016_baroclinic_wave
    CASE ('dcmip2016_cyclone')
        CALL getin_etat0_dcmip2016_cyclone
    CASE ('dcmip2016_supercell')
        CALL getin_etat0_dcmip2016_supercell
    CASE DEFAULT
       collocated=.FALSE.
       autoinit_mass = .FALSE.
    END SELECT

    !------------------- Ad hoc interfaces --------------------
    SELECT CASE (TRIM(etat0_type))
     CASE ('database')
        CALL etat0_database(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
    CASE ('start_file')
       CALL etat0_start_file(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 ('held_suarez')
       PRINT *,"Held & Suarez (1994) test case"
       CALL etat0_heldsz(f_ps,f_phis,f_theta_rhodz,f_u, f_q)
    CASE ('venus')
       CALL etat0_venus(f_ps, f_phis, f_theta_rhodz, f_u, f_q)
       PRINT *, "Venus (Lebonnois et al., 2012) test case"
   CASE DEFAULT
      IF(collocated) THEN
         CALL etat0_collocated(f_phis,f_ps,f_mass,f_theta_rhodz,f_u, f_geopot,f_W, f_q)
      ELSE
         PRINT*, 'Bad selector for variable etat0 <',TRIM(etat0_type),'>'// &
            ' options are  <isothermal>, <temperature_profile>, <jablonowsky06>, <dcmip5>, <williamson91.6>,'& 
                         //' <start_file>, <database>, <academic>, <held_suarez>, <venus>, <dcmip1>,'         &
                         //' <dcmip2_mountain,dcmip2_schaer_noshear,dcmip2_schaer_shear>, <dcmip3>, <dcmip4>,'&
                         //' <dcmip2016_baroclinic_wave>, <dcmip2016_cyclone>, <dcmip2016_supercell>'
         STOP
      END IF
    END SELECT

    IF(autoinit_mass) THEN
       DO ind=1,ndomain
          IF (.NOT. assigned_domain(ind) .OR. .NOT. is_omp_level_master) CYCLE
          CALL swap_dimensions(ind)
          CALL swap_geometry(ind)
          mass=f_mass(ind); ps=f_ps(ind)
          CALL compute_rhodz(.TRUE., ps, mass) ! initialize mass distribution using ps
       END DO
    END IF
 
  END SUBROUTINE etat0

END MODULE etat0_mod