MODULE etat0_mod
   USE icosa
   USE etat0_collocated_mod
   IMPLICIT NONE
   PRIVATE

   REAL(rstd) :: etat0_temp

   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

      ! @getin_doc etat0 realm=initial type=String default=none
      !  Selects initial condition for dynamics : surface geopotential, surface pressure, wind,
      ! temperature, composition. Valid values for idealized test cases:
      !  * 'williamson91.6' : Williamson (1991) test case 6 (rotating shallow-water Rossby-Haurwitz wave)
      !  * `isothermal` : isothermal atmosphere at rest, flat ground.
      !  * `temperature profile` : user-profided vertical temperature profile, rest, flat.
      !  * `jablonowsky06` : Jablonowsky & Williamson (2006) baroclinic instabilty.
      !  * `dcmip1` : DCMIP 2012 test case 1
      !  * `dcmip2_mountain` : DCMIP 2012 test case 2 (state of rest with orography)
      !  * `dcmip2_schaer_noshear` : DCMIP 2012 test case 2 (orographic wave with non-sheared wind)
      !  * `dcmip2_schaer_shear` : DCMIP 2012 test case 2 (orographic wave with sheared wind)
      !  * `dcmip3` : DCMIP 2012 test case 3
      !  * `dcmip4` : DCMIP 2012 test case 4 (baroclinic instability)
      !  * `dcmip5` : DCMIP 2012 test case 5 (cyclone)
      !  * `venus`  : Venus-like
      !  * 'dcmip2016_baroclinic_wave'
      !  * 'dcmip2016_cyclone'
      !  * 'dcmip2016_supercell'
      !  * 'academic'
      !  * 'held_suarez'
      !  * 'bubble' : zonally-symmetric warm bubble
      ! Valid values for realistic initial conditions, read from external file(s) via XIOS :
      !  * 'start_file'
      !  * 'start_file_interp'
      !  * 'database'

      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
