MODULE etat0_dcmip1_mod
   USE icosa
   IMPLICIT NONE
   PRIVATE

   INTEGER, PARAMETER :: const = 1, cos_bell = 2, slotted_cyl = 3, &
                         dbl_cos_bell_q1 = 4, dbl_cos_bell_q2 = 5, complement = 6, hadley = 7, &
                         dcmip11 = -1
   INTEGER, SAVE :: shape
!$OMP THREADPRIVATE(shape)

   REAL(rstd), SAVE  :: h0 = 1.
!$OMP THREADPRIVATE(h0)
   REAL(rstd), SAVE  :: lon0 = 3*pi/2
!$OMP THREADPRIVATE(lon0)
   REAL(rstd), SAVE  :: lat0 = 0.0
!$OMP THREADPRIVATE(lat0)
   REAL(rstd), SAVE  :: R0
!$OMP THREADPRIVATE(R0)
   REAL(rstd), SAVE  :: latc1 = 0.
!$OMP THREADPRIVATE(latc1)
   REAL(rstd), SAVE  :: latc2 = 0.
!$OMP THREADPRIVATE(latc2)
   REAL(rstd), SAVE  :: lonc1 = 5*pi/6
!$OMP THREADPRIVATE(lonc1)
   REAL(rstd), SAVE  :: lonc2 = 7*pi/6
!$OMP THREADPRIVATE(lonc2)
   REAL(rstd), SAVE  :: zt = 1000.0
!$OMP THREADPRIVATE(zt)
   REAL(rstd), SAVE  :: rt
!$OMP THREADPRIVATE(rt)
   REAL(rstd), SAVE  :: zc = 5000.0
!$OMP THREADPRIVATE(zc)

   PUBLIC getin_etat0, compute_etat0

CONTAINS

   SUBROUTINE getin_etat0
      CHARACTER(len=255) :: dcmip1_adv_shape
      R0 = radius*0.5
      rt = radius*0.5
      dcmip1_adv_shape = 'cos_bell'
      CALL getin('dcmip1_shape', dcmip1_adv_shape)
      SELECT CASE (TRIM(dcmip1_adv_shape))
      CASE ('const')
         shape = const
      CASE ('cos_bell')
         shape = cos_bell
      CASE ('slotted_cyl')
         shape = slotted_cyl
      CASE ('dbl_cos_bell_q1')
         shape = dbl_cos_bell_q1
      CASE ('dbl_cos_bell_q2')
         shape = dbl_cos_bell_q2
      CASE ('complement')
         shape = complement
      CASE ('hadley')  ! hadley like meridional circulation
         shape = hadley
      CASE ('dcmip11')
         IF (nqtot < 5) THEN
            PRINT *, 'Error : etat0_dcmip=dcmip11 and nqtot = ', nqtot, ' .'
            PRINT *, 'nqtot must be equal to 5 when etat0_dcmip=dcmip11'
            STOP
         END IF
         shape = dcmip11
      CASE DEFAULT
         PRINT *, 'Bad selector for variable dcmip1_adv_shape : <', TRIM(dcmip1_adv_shape), &
            '> options are <const>, <slotted_cyl>, <cos_bell>, <dbl_cos_bell_q1>', &
            '<dbl_cos_bell_q2>, <complement>, <hadley>'
         STOP
      END SELECT
   END SUBROUTINE getin_etat0

   SUBROUTINE compute_etat0(ngrid, lon, lat, phis, ps, temp, ulon, ulat, q)
      INTEGER, INTENT(IN) :: ngrid
      REAL(rstd), INTENT(IN) :: lon(ngrid), lat(ngrid)
      REAL(rstd), INTENT(OUT) :: ps(ngrid), phis(ngrid)
      REAL(rstd), INTENT(OUT) :: temp(ngrid, llm), ulon(ngrid, llm), ulat(ngrid, llm)
      REAL(rstd), INTENT(OUT) :: q(ngrid, llm, nqtot)
      ps = ncar_p0
      phis = 0.; temp = 0.; 
      ulon = 0.; ulat = 0.
      SELECT CASE (shape)
      CASE (dcmip11)
         CALL compute_etat0_ncar(4, ngrid, lon, lat, q(:, :, 1))
         CALL compute_etat0_ncar(5, ngrid, lon, lat, q(:, :, 2))
         CALL compute_etat0_ncar(3, ngrid, lon, lat, q(:, :, 3))
         CALL compute_etat0_ncar(6, ngrid, lon, lat, q(:, :, 4))
         CALL compute_etat0_ncar(1, ngrid, lon, lat, q(:, :, 5))
      CASE DEFAULT
         CALL compute_etat0_ncar(shape, ngrid, lon, lat, q(:, :, 1))
      END SELECT
   END SUBROUTINE compute_etat0

   SUBROUTINE compute_etat0_ncar(icase, ngrid, lon, lat, q)
      USE disvert_mod
      USE omp_para
      INTEGER, INTENT(IN) :: icase, ngrid
      REAL(rstd), INTENT(IN) :: lon(ngrid), lat(ngrid)
      REAL(rstd), INTENT(OUT) :: q(ngrid, llm)
      REAL(rstd) :: zr(llm + 1), zrl(llm), qxt1(ngrid, llm)
      REAL(rstd) :: pr
      !  REAL(rstd) :: lon, lat
      INTEGER :: l

      DO l = 1, llm + 1
         pr = ap(l) + bp(l)*ncar_p0
         zr(l) = -kappa*cpp*ncar_T0/g*log(pr/ncar_p0)
      END DO

      DO l = 1, llm
         zrl(l) = 0.5*(zr(l) + zr(l + 1))
      END DO

      SELECT CASE (icase)
      CASE (1)
         q = 1
      CASE (2)
         CALL cosine_bell_1(q)
      CASE (3)
         CALL slotted_cylinders(q)
      CASE (4)
         CALL cosine_bell_2(q)
      CASE (5)
         CALL cosine_bell_2(q)
         DO l = 1, llm
            q(:, l) = 0.9 - 0.8*q(:, l)*q(:, l)
         END DO
      CASE (6)
         ! tracer such that, in combination with the other tracer fields
         ! with weight (3/10), the sum is equal to one
         CALL cosine_bell_2(qxt1)
         DO l = 1, llm
            q(:, l) = 0.9 - 0.8*qxt1(:, l)*qxt1(:, l)
         END DO
         q = q + qxt1
         CALL slotted_cylinders(qxt1)
         q = q + qxt1
         q = 1.-q*0.3
      CASE (7)  ! hadley like meridional circulation
         CALL hadleyq(q)
      END SELECT

   CONTAINS
!======================================================================

      SUBROUTINE cosine_bell_1(hx)
         REAL(rstd) :: hx(ngrid, llm)
         REAL(rstd) :: rr1
         INTEGER :: n, l
         DO l = ll_begin, ll_end
            DO n = 1, ngrid
               CALL dist_lonlat(lon0, lat0, lon(n), lat(n), rr1)  ! GC distance from center
               rr1 = radius*rr1
               IF (rr1 .LT. R0) then
                  hx(n, l) = 0.5*h0*(1 + cos(pi*rr1/R0))
               ELSE
                  hx(n, l) = 0.0
               END IF
            END DO
         END DO
      END SUBROUTINE cosine_bell_1

!==============================================================================

      SUBROUTINE cosine_bell_2(hx)
         REAL(rstd) :: hx(ngrid, llm)
         REAL(rstd) :: rr1, rr2, dd1, dd2, dd1t1, dd1t2, dd2t1
         INTEGER :: n, l
         DO l = ll_begin, ll_end
            DO n = 1, ngrid
               CALL dist_lonlat(lonc1, latc1, lon(n), lat(n), rr1)  ! GC distance from center
               rr1 = radius*rr1
               CALL dist_lonlat(lonc2, latc2, lon(n), lat(n), rr2)  ! GC distance from center
               rr2 = radius*rr2
               dd1t1 = rr1/rt
               dd1t1 = dd1t1*dd1t1
               dd1t2 = (zrl(l) - zc)/zt
               dd1t2 = dd1t2*dd1t2
               dd1 = dd1t1 + dd1t2
               dd1 = Min(1.0, dd1)
               dd2t1 = rr2/rt
               dd2t1 = dd2t1*dd2t1
               dd2 = dd2t1 + dd1t2
               dd2 = Min(1.0, dd2)
               hx(n, l) = 0.5*(1.+cos(pi*dd1)) + 0.5*(1.+cos(pi*dd2))
            END DO
         END DO
      END SUBROUTINE cosine_bell_2

      !=============================================================================

      SUBROUTINE slotted_cylinders(hx)
         REAL(rstd) :: hx(ngrid, llm)
         REAL(rstd) :: rr1, rr2, dd1, dd2, dd1t1, dd1t2, dd2t1
         INTEGER :: n, l
         DO l = ll_begin, ll_end
            DO n = 1, ngrid
               CALL dist_lonlat(lonc1, latc1, lon(n), lat(n), rr1)  ! GC distance from center
               rr1 = radius*rr1
               CALL dist_lonlat(lonc2, latc2, lon(n), lat(n), rr2)  ! GC distance from center
               rr2 = radius*rr2
               dd1t1 = rr1/rt
               dd1t1 = dd1t1*dd1t1
               dd1t2 = (zrl(l) - zc)/zt
               dd1t2 = dd1t2*dd1t2
               dd1 = dd1t1 + dd1t2
               dd2t1 = rr2/rt
               dd2t1 = dd2t1*dd2t1
               dd2 = dd2t1 + dd1t2
               IF (dd1 .LT. 0.5) Then
                  hx(n, l) = 1.0
               ELSEIF (dd2 .LT. 0.5) Then
                  hx(n, l) = 1.0
               ELSE
                  hx(n, l) = 0.1
               END IF
               IF (zrl(l) .GT. zc) Then
                  IF (ABS(latc1 - lat_i(n)) .LT. 0.125) Then
                     hx(n, l) = 0.1
                  END IF
               END IF
            END DO
         END DO
      END SUBROUTINE slotted_cylinders

!==============================================================================

      SUBROUTINE hadleyq(hx)
         REAL(rstd)::hx(ngrid, llm)
         REAL(rstd), PARAMETER:: zz1 = 2000., zz2 = 5000., zz0 = 0.5*(zz1 + zz2)
         INTEGER :: l

         DO l = ll_begin, ll_end
            IF ((zz1 .LT. zrl(l)) .and. (zrl(l) .LT. zz2)) THEN
               hx(:, l) = 0.5*(1.+cos(2*pi*(zrl(l) - zz0)/(zz2 - zz1)))
            ELSE
               hx(:, l) = 0.0
            END IF
         END DO
      END SUBROUTINE hadleyq

   END SUBROUTINE compute_etat0_ncar

END MODULE etat0_dcmip1_mod
