MODULE guided_ncar_mod
  USE icosa
  PRIVATE
  
  INTEGER,SAVE :: case_wind
!$OMP THREADPRIVATE(case_wind)
  
  REAL(rstd), PARAMETER :: alpha=0.0 ! tilt of solid-body rotation
  REAL(rstd), PARAMETER :: tau_hadley=daysec, tau = 12*daysec ! 12 days               ! see p. 16
  REAL(rstd), PARAMETER :: w0_deform = 23000*pi/tau, b=0.2, ptop=25494.4  ! see p. 16 
  REAL(rstd), PARAMETER :: u0_hadley=40.,w0_hadley=0.15 ,ztop= 12000. 
  INTEGER, PARAMETER    :: K_hadley=5

  PUBLIC init_guided, guided

CONTAINS

  SUBROUTINE init_guided
    IMPLICIT NONE
    CHARACTER(LEN=255) :: wind
    wind='deform'
    CALL getin('dcmip1_wind',wind)
    SELECT CASE(TRIM(wind))
    CASE('solid')
       case_wind=0
    CASE('deform')
       case_wind=1
    CASE('hadley')
       case_wind=2
    CASE DEFAULT
       PRINT*,'Bad selector for variable ncar_adv_wind : <', TRIM(wind),'> options are <solid>, <deform>, <hadley>'
       END SELECT
  END SUBROUTINE init_guided
  
  SUBROUTINE guided(tt, f_ps, f_theta_rhodz, f_u, f_q)
  USE icosa
    IMPLICIT NONE
    REAL(rstd), INTENT(IN):: tt
    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 :: ue(:,:)
    INTEGER :: ind
    
    DO ind = 1 , ndomain
      IF (.NOT. assigned_domain(ind)) CYCLE
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind) 
      ue = f_u(ind)  
      CALL wind_profile(tt,ue)
    END DO    

  END SUBROUTINE guided
  

  SUBROUTINE wind_profile(tt,ue)
    USE icosa
    USE disvert_mod
    IMPLICIT NONE
    REAL(rstd),INTENT(IN)  :: tt ! current time 
    REAL(rstd),INTENT(OUT) :: ue(iim*3*jjm,llm)
    REAL(rstd) :: lon, lat
    REAL(rstd) :: nx(3),n_norm,Velocity(3,llm)
    REAL(rstd) :: rr1,rr2,bb,cc,aa,hmx
    REAL(rstd) :: v1(3),v2(3),ny(3)
    INTEGER :: i,j,n,l
    REAL(rstd) :: pitbytau,kk, pr, zr, u0, u1, v0

    pitbytau = pi*tt/tau
    kk = 10*radius/tau
    u0 = 2*pi*radius/tau ! for solid-body rotation
!---------------------------------------------------------
    DO l = 1,llm 
       pr = presnivs(l)
       zr = -kappa*cpp*ncar_T0/g*log(pr/ncar_p0)  ! reciprocal of (1) p. 13, isothermal atmosphere
       u1 = w0_deform*radius/b/ptop*cos(2*pitbytau)*(exp((ptop-pr)/b/ptop)-exp((pr-ncar_p0)/b/ptop))
       v0 = -radius*w0_hadley*pi/(5.0*ztop)*(ncar_p0/pr)*cos(pi*zr/ztop)*cos(pi*tt/tau_hadley) ! for Hadley cell

       DO j=jj_begin-1,jj_end+1
          DO i=ii_begin-1,ii_end+1
             n=(j-1)*iim+i
             CALL compute_velocity(xyz_e(n+u_right,:),l,velocity(:,l))
             CALL cross_product2(xyz_v(n+z_rdown,:)/radius,xyz_v(n+z_rup,:)/radius,nx)       
             ue(n+u_right,l)=1e-10
             n_norm=sqrt(sum(nx(:)**2))
             IF (n_norm>1e-30) THEN
                nx=-nx/n_norm*ne(n,right)
                ue(n+u_right,l)=sum(nx(:)*velocity(:,l))
                IF (ABS(ue(n+u_right,l))<1e-100) PRINT *,"ue(n+u_right)==0",i,j,velocity(:,1)
             ENDIF
             
             CALL compute_velocity(xyz_e(n+u_lup,:),l,velocity(:,l))
             CALL cross_product2(xyz_v(n+z_up,:)/radius,xyz_v(n+z_lup,:)/radius,nx)
             
             ue(n+u_lup,l)=1e-10
             n_norm=sqrt(sum(nx(:)**2))
             IF (n_norm>1e-30) THEN
                nx=-nx/n_norm*ne(n,lup)
                ue(n+u_lup,l)=sum(nx(:)*velocity(:,l))
             ENDIF
             
             CALL compute_velocity(xyz_e(n+u_ldown,:),l,velocity(:,l))
             CALL cross_product2(xyz_v(n+z_ldown,:)/radius,xyz_v(n+z_down,:)/radius,nx)
             
             ue(n+u_ldown,l)=1e-10
             n_norm=sqrt(sum(nx(:)**2))
             IF (n_norm>1e-30) THEN
                nx=-nx/n_norm*ne(n,ldown)
                ue(n+u_ldown,l)=sum(nx(:)*velocity(:,l))
                IF (ABS(ue(n+u_ldown,l))<1e-100) PRINT *,"ue(n+u_ldown)==0",i,j
             ENDIF
          ENDDO
       ENDDO
    END DO
    
  CONTAINS
	   
    SUBROUTINE compute_velocity(x,l,velocity)
      IMPLICIT NONE
      REAL(rstd),INTENT(IN)  :: x(3)
      INTEGER,INTENT(IN)::l
      REAL(rstd),INTENT(OUT) :: velocity(3)
      REAL(rstd) :: e_lat(3), e_lon(3)
      REAL(rstd) :: lon,lat
      REAL(rstd) :: u,v
	      
      CALL xyz2lonlat(x/radius,lon,lat)
      e_lat(1) = -cos(lon)*sin(lat)
      e_lat(2) = -sin(lon)*sin(lat)
      e_lat(3) = cos(lat)
        
      e_lon(1) = -sin(lon)
      e_lon(2) = cos(lon)
      e_lon(3) = 0

      u = 0.0 ; v = 0.0
    
      SELECT CASE(case_wind) 
      CASE(0)  ! Solid-body rotation
         u=u0*(cos(lat)*cos(alpha)+sin(lat)*sin(alpha)*cos(lon))
         v=-u0*sin(lon)*sin(alpha)
      CASE(1)  ! 3D Deformational flow - 
         lon = lon-2*pitbytau
         u = kk*sin(lon)*sin(lon)*sin(2*lat)*cos(pitbytau)+ u0*cos(lat)
         u = u + u1*cos(lon)*cos(lat)**2
         v = kk*sin(2*lon)*cos(lat)*cos(pitbytau) 
      CASE(2) ! Hadley-like flow
         u = u0_hadley*cos(lat)
         v = v0*cos(lat)*sin(5.*lat)    ! Eq. 37 p. 19
      CASE DEFAULT 
         PRINT*,"not valid choice of wind" 
      END SELECT
      
      Velocity(:)=(u*e_lon(:)+v*e_lat(:)+1e-50)
      
    END SUBROUTINE compute_velocity
    
  END SUBROUTINE  wind_profile


END MODULE guided_ncar_mod