source: codes/icosagcm/trunk/src/guided_ncar_mod.f90 @ 19

MODULE guided_ncar_mod
  USE icosa
  PRIVATE
  
  REAL(rstd),SAVE :: dt
  INTEGER,SAVE    :: test
  
  PUBLIC init_guided, guided

CONTAINS


  SUBROUTINE init_guided(dt_in)
  IMPLICIT NONE
    REAL(rstd),INTENT(IN) :: dt_in
    
    dt=dt_in
    test=2
    
  END SUBROUTINE init_guided

  
  SUBROUTINE guided(it, f_ps, f_theta_rhodz, f_u, f_q)
  USE icosa
    IMPLICIT NONE
    INTEGER, INTENT(IN)   :: it
    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
      CALL swap_dimensions(ind)
      CALL swap_geometry(ind) 
      ue = f_u(ind)  
      CALL wind_profile(it,ue)
    END DO    

  END SUBROUTINE guided
  

  SUBROUTINE wind_profile(it,ue)
  USE icosa
  IMPLICIT NONE
    INTEGER,INTENT(IN) :: it 
    REAL(rstd),INTENT(OUT) :: ue(iim*3*jjm,llm)
    REAL(rstd) :: lon, lat
    REAL(rstd) :: nx(3),n_norm,Velocity(3,llm)
    REAL(rstd), PARAMETER :: lon0=3*pi/2,lat0=0.0
    REAL(rstd) :: rr1,rr2,bb,cc,aa,hmx
    REAL(rstd) :: v1(3),v2(3),ny(3)
    INTEGER :: i,j,n,l
    REAL(rstd) :: zrl(llm)

    CALL compute_zrl(zrl)
!---------------------------------------------------------
    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 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) :: u0
      REAL(rstd)  :: e_lat(3), e_lon(3)
      REAL(rstd) :: lon,lat
      REAL(rstd),PARAMETER::alpha=0.0
      REAL(rstd) :: u,v
!--------------------------------- test 1
      REAL(rstd),parameter:: nday=12,daysec=86400
      REAL(rstd),PARAMETER :: tau = nday*daysec
      REAL(rstd) :: tt, pitbytau,kk
!---------------------------------- hadley 
      INTEGER,PARAMETER::K_hadly=5
      REAL(rstd),PARAMETER::u0_hadly=40.,w0_hadly=0.25 ,ztop= 12000. 
      REAL(rstd)::coff_hadly
        
        
      tt = it*dt 
      pitbytau = pi*tt/tau 
      kk = 10*radius/tau 
      
      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
      u0 = 2*pi*radius/tau
    
      SELECT CASE(test) 
        CASE(0)  
          u=u0*(cos(lat)*cos(alpha)+sin(lat)*sin(alpha)*cos(lon))
          v=-u0*sin(lon)*sin(alpha)
        CASE(1) 
          lon = lon - 2*pitbytau 
          u= kk*sin(lon)*sin(lon)*sin(2*lat)*cos(pitbytau)+ u0*cos(lat) 
          v= kk*sin(2*lon)*cos(lat)*cos(pitbytau) 
        CASE(2)
          coff_hadly = -radius*w0_hadly*pi/(5.0*ztop)  
          u= u0_hadly*cos(lat) 
          v= coff_hadly*cos(lat)*sin(5.*lat)*cos(pi*zrl(l)/ztop)*cos(pitbytau)
        CASE DEFAULT 
          PRINT*,"not valid choice of wind" 
       END SELECT 

      Velocity(:)=(u*e_lon(:)+v*e_lat(:)+1e-50)

    END SUBROUTINE compute_velocity
    
    SUBROUTINE compute_zrl(zrl)
    USE disvert_mod
    IMPLICIT NONE
      REAL(rstd),INTENT(OUT) :: zrl(llm)
      INTEGER :: l
      REAL(rstd) :: zr(llm+1)
      REAl(rstd) :: pr
      REAL(rstd), PARAMETER :: T0=300
      REAL(rstd), PARAMETER :: R_d=287.0   
      REAL(rstd), PARAMETER :: ps0=100000.0

      
      DO    l    = 1, llm+1
        pr = ap(l) + bp(l)*ps0
        zr(l) = -R_d*T0/g*log(pr/ps0) 
      ENDDO    

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


END MODULE guided_ncar_mod