MODULE etat0_ncar_mod
  USE icosa
  PRIVATE

  REAL(rstd), SAVE  :: h0=1.
  REAL(rstd), SAVE  :: lon0=3*pi/2
  REAL(rstd), SAVE  :: lat0=0.0
  REAL(rstd), SAVE  :: alpha=0.0
  REAL(rstd), SAVE  :: R0 
  REAL(rstd), SAVE  :: lat1=0.
  REAL(rstd), SAVE  :: lat2=0.
  REAL(rstd), SAVE  :: lon1=pi/6
  REAL(rstd), SAVE  :: lon2=-pi/6
  REAL(rstd), SAVE  :: latc1=0.
  REAL(rstd), SAVE  :: latc2=0.
  REAL(rstd), SAVE  :: lonc1=5*pi/6
  REAL(rstd), SAVE  :: lonc2=7*pi/6
  REAL(rstd), SAVE  :: zt=1000.0
  REAL(rstd), SAVE  :: rt
  REAL(rstd), SAVE  :: zc=5000.0

  PUBLIC etat0
  
CONTAINS
  
  SUBROUTINE etat0(f_ps,f_phis,f_theta_rhodz,f_u,f_q)
  USE icosa
  IMPLICIT NONE
    TYPE(t_field),POINTER :: f_ps(:)
    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(:)
    REAL(rstd),POINTER :: phis(:)
    REAL(rstd),POINTER :: theta_rhodz(:,:)
    REAL(rstd),POINTER :: u(:,:)
    REAL(rstd),POINTER :: q(:,:,:)
    INTEGER :: ind
    
    R0=radius*0.5
    rt=radius*0.5
    
    DO ind=1,ndomain
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind)
      ps=f_ps(ind)
      phis=f_phis(ind)
      theta_rhodz=f_theta_rhodz(ind)
      u=f_u(ind)
      q=f_q(ind)
      CALL compute_etat0_ncar(ps, phis, theta_rhodz, u, q(:,:,1))
    ENDDO

  END SUBROUTINE etat0
  
  SUBROUTINE compute_etat0_ncar(ps, phis, theta_rhodz, u, q)
  USE icosa
  USE disvert_mod
  USE pression_mod
  USE exner_mod
  USE geopotential_mod
  USE theta2theta_rhodz_mod
  IMPLICIT NONE  
  REAL(rstd),INTENT(OUT) :: ps(iim*jjm)
  REAL(rstd),INTENT(OUT) :: phis(iim*jjm)
  REAL(rstd),INTENT(OUT) :: theta_rhodz(iim*jjm,llm)
  REAL(rstd),INTENT(OUT) :: u(3*iim*jjm,llm)
  REAL(rstd),INTENT(OUT) :: q(iim*jjm,llm)  


  REAL(rstd) :: qxt1(iim*jjm,llm)
  REAL(rstd) :: lon, lat
  REAL(rstd) ::dd1,dd2,dd1t1,dd1t2,dd2t1
  REAL(rstd) :: pr, zr(llm+1), zrl(llm)
  REAL(rstd) :: rr1,rr2,bb,cc,aa,hmx
  REAL(rstd) :: X2(3),X1(3)
  INTEGER :: i,j,n,l
  CHARACTER(len=255) :: ncar_adv_shape

  u = 0.0 ; phis = 0 ; theta_rhodz = 0 ; ps = ncar_p0
  
  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
  
  ncar_adv_shape='cos_bell'
  CALL getin('ncar_adv_shape',ncar_adv_shape)
  
  SELECT CASE(TRIM(ncar_adv_shape))
     !--------------------------------------------- SINGLE COSINE BELL
  CASE('const')
     q=1
  CASE('cos_bell')
     CALL cosine_bell_1(q) 
     
  CASE('slotted_cyl')
     CALL slotted_cylinders(q) 	
     
  CASE('dbl_cos_bell_q1')
     CALL cosine_bell_2(q) 
     
  CASE('dbl_cos_bell_q2')
     CALL cosine_bell_2(q) 
     DO l=1,llm
        q(:,l)= 0.9 - 0.8*q(:,l)*q(:,l)
     END DO
     
  CASE('complement')
     ! 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('hadley')  ! hadley like meridional circulation 
     CALL hadleyq(q) 
     
  CASE DEFAULT
     PRINT *, 'Bad selector for variable ncar_adv_shape : <', TRIM(ncar_adv_shape),   &
          '> options are <const>, <slotted_cyl>, <cos_bell>, <dbl_cos_bell_q1>', &
          '<dbl_cos_bell_q2>, <complement>, <hadley>'
     STOP
     
  END SELECT
      
CONTAINS

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

    SUBROUTINE cosine_bell_1(hx)
    IMPLICIT NONE 
    REAL(rstd) :: hx(iim*jjm,llm)
    
      DO l=1,llm 
        DO j=jj_begin-1,jj_end+1
          DO i=ii_begin-1,ii_end+1
            n=(j-1)*iim+i
            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
            CALL dist_lonlat(lon0,lat0,lon,lat,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 DO
   
    END SUBROUTINE cosine_bell_1

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

    SUBROUTINE   cosine_bell_2(hx) 
    IMPLICIT NONE 	
    REAL(rstd) :: hx(iim*jjm,llm)
   
      DO l=1,llm
        DO j=jj_begin-1,jj_end+1
          DO i=ii_begin-1,ii_end+1
            n=(j-1)*iim+i
            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center 
            rr1 = radius*rr1
            CALL dist_lonlat(lonc2,latc2,lon,lat,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 DO

    END SUBROUTINE cosine_bell_2

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

    SUBROUTINE slotted_cylinders(hx) 
    IMPLICIT NONE
    REAL(rstd) :: hx(iim*jjm,llm)
    
      DO l=1,llm
        DO j=jj_begin-1,jj_end+1
          DO i=ii_begin-1,ii_end+1
            n=(j-1)*iim+i
            CALL xyz2lonlat(xyz_i(n,:),lon,lat)
            CALL dist_lonlat(lonc1,latc1,lon,lat,rr1)  ! GC distance from center 
            rr1 = radius*rr1
            CALL dist_lonlat(lonc2,latc2,lon,lat,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) .LT. 0.125 ) Then 
                hx(n,l)= 0.0
 	      ENDIF 
            ENDIF  
  
          ENDDO
       END DO
      END DO

    END SUBROUTINE slotted_cylinders 

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

    SUBROUTINE hadleyq(hx) 
    IMPLICIT NONE 
    REAL(rstd)::hx(iim*jjm,llm) 
    REAL(rstd),PARAMETER:: zz1=3500.,zz2=6500.,zz0=0.5*(zz1+zz2)
   
      DO l=1,llm
        IF ( ( zz1 .LT. zrl(l) ) .and. ( zrl(l) .LT. zz2 ) )  THEN 
          hx(:,l) = 0.5*(1. + cos(0.002*pi*(zrl(l)-zz0)/3.))  
        ELSE
	  hx(:,l) = 0.0 
        END IF
      END DO
    END SUBROUTINE hadleyq 

  END SUBROUTINE compute_etat0_ncar
  

END MODULE etat0_ncar_mod