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  :: alpha=0.0
!$OMP THREADPRIVATE(alpha)
  REAL(rstd), SAVE  :: R0 
!$OMP THREADPRIVATE(R0)
  REAL(rstd), SAVE  :: lat1=0.
!$OMP THREADPRIVATE(lat1)
  REAL(rstd), SAVE  :: lat2=0.
!$OMP THREADPRIVATE(lat2)
  REAL(rstd), SAVE  :: lon1=pi/6
!$OMP THREADPRIVATE(lon1)
  REAL(rstd), SAVE  :: lon2=-pi/6
!$OMP THREADPRIVATE(lon2)
  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 :: n,l
  
  DO l=1, llm+1
     pr = ap(l) + bp(l)*ncar_p0
     zr(l) = -kappa*cpp*ncar_T0/g*log(pr/ncar_p0) 
  ENDDO
  
  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,rr2   
    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
             ENDIF
          ENDIF
       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 :: n,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