meteo.f90

Go to the documentation of this file.

!********************************************************************
!
!********************************************************************

MODULE meteo_module

  USE variables, ONLY: gi   ! Grav acceleration 

  IMPLICIT NONE

  REAL*8 :: gt     
  
  REAL*8 :: gs    
  
  REAL*8 :: p0     
  
  REAL*8 :: t0     
  
  REAL*8 :: h1    
  
  REAL*8 :: h2    
  
  REAL*8 :: u_atm0
  
  REAL*8 :: duatm_dz0  

  REAL*8 :: cos_theta , sin_theta

  REAL*8 :: rho_atm0

  REAL*8 :: u_atm   

  REAL*8 :: rho_atm  

  REAL*8 :: visc_atm 

  REAL*8 :: visc_atm0

  REAL*8 :: ta      

  REAL*8 :: pa      

  REAL*8 :: dpdz     

  REAL*8 :: dtdz     

  REAL*8 :: duatm_dz 

  REAL :: rair

  REAL :: cpair

  INTEGER :: n_atm_profile

  REAL*8, ALLOCATABLE :: atm_profile(:,:)

  CHARACTER*10 :: read_atm_profile

  REAL*8 :: u_r , z_r , exp_wind

  REAL*8, ALLOCATABLE :: rho_atm_month_lat(:) , pres_atm_month_lat(:) ,   &
       temp_atm_month_lat(:) , temp_atm_month(:,:)

  REAL*8, ALLOCATABLE :: h_levels(:)
  
  REAL*8 :: wind_mult_coeff

  SAVE

CONTAINS

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE initialize_meteo

    IMPLICIT NONE

    CALL zmet

    rho_atm0 = rho_atm

    RETURN

  END SUBROUTINE initialize_meteo

  !******************************************************************************
  !
  !******************************************************************************

  SUBROUTINE zmet

    USE plume_module, ONLY: z , vent_height

    USE variables, ONLY : verbose_level

    IMPLICIT NONE

    REAL*8 :: const, const1, t1, p1, p2
    REAL*8 :: const2

    REAL*8 :: we_wind , ns_wind

    REAL*8 :: cs

    IF ( read_atm_profile .EQ. 'card' ) THEN

       ! interp density profile
       CALL interp_1d_scalar(atm_profile(1,:), atm_profile(2,:), z, rho_atm)

       ! interp pressure profile
       CALL interp_1d_scalar(atm_profile(1,:), atm_profile(3,:), z, pa)

       ! interp temperature profile
       CALL interp_1d_scalar(atm_profile(1,:), atm_profile(4,:), z, ta)

       ! interp pressure profile
       CALL interp_1d_scalar(atm_profile(1,:), atm_profile(6,:), z, we_wind)

       ! interp pressure profile
       CALL interp_1d_scalar(atm_profile(1,:), atm_profile(7,:), z, ns_wind)

       u_atm = dsqrt( we_wind**2 + ns_wind**2 )

       cos_theta = we_wind / u_atm
       sin_theta = ns_wind / u_atm

    ELSEIF ( read_atm_profile .EQ. 'table' ) THEN

       ! interp density profile
       CALL interp_1d_scalar(h_levels(:), rho_atm_month_lat(:), z, rho_atm)

       ! interp pressure profile
       CALL interp_1d_scalar(h_levels(:), pres_atm_month_lat(:), z, pa)

       ! interp temperature profile
       CALL interp_1d_scalar(h_levels(:), temp_atm_month_lat(:), z, ta)

       IF ( z .LE. z_r ) THEN

          u_atm = u_r * ( ( z - vent_height ) / z_r )**exp_wind

       ELSE

          u_atm = u_r

       END IF

       cos_theta = 1.d0
       sin_theta = 0.d0


    ELSEIF ( read_atm_profile .EQ. 'standard' ) THEN


       ! u_atm = u_atm0 + duatm_dz * z 
       ! duatm_dz = duatm_dz0
       
       IF ( z .LE. z_r ) THEN

          u_atm = u_r * ( z / z_r )**exp_wind

       ELSE

          u_atm = u_r

       END IF

       cos_theta = 1.d0
       sin_theta = 0.d0


       !      
       ! ... Temperature and pressure at the tropopause bottom
       !
       const = - gi / ( rair * gt )
       t1 = t0 + gt * h1
       p1 = p0 * (t1/t0)**const
       const1 = gi / ( rair * t1 )
       p2 = p1 * dexp( -const1 * ( h2-h1 ) )
       const2 = - gi / ( rair * gs )

       IF ( z <= h1 ) THEN

          ! ... Troposphere

          ta = t0 + gt * z
          pa = p0 * ( ta / t0 )**const
          rho_atm = pa / ( rair*ta )

       ELSE IF (z > h1 .AND. z <= h2) THEN

          ! ... Tropopause

          ta = t1
          pa = p1 * dexp( -const1 * ( z - h1 ) )
          rho_atm = pa / ( rair * ta )

       ELSE

          ! ... Stratosphere

          ta = t0 + gt*h1 + gs*(z-h2)
          pa = p2 * (ta/t1)**const2 
          rho_atm = pa / (rair*ta)

       ENDIF

    END IF

    ! ... Air viscosity ( Armienti et al. 1988)
    cs = 120.d0
    visc_atm = visc_atm0 * ( 288.15d0 + cs ) / ( ta + cs ) * ( ta / 288.15d0 )**1.5d0

    IF ( verbose_level .GE. 1 ) THEN

       WRITE(*,*) z,rho_atm,pa,ta,u_atm,cos_theta,sin_theta
       WRITE(*,*) 'visc_atm',visc_atm
       READ(*,*)

    END IF

    RETURN

  END SUBROUTINE zmet

!---------------------------------------------------------------------------
!
!---------------------------------------------------------------------------

  SUBROUTINE interp_1d_scalar(x1, f1, x2, f2)
    IMPLICIT NONE
    
    REAL*8, INTENT(IN), DIMENSION(:) :: x1, f1
    REAL*8, INTENT(IN) :: x2
    REAL*8, INTENT(OUT) :: f2
    INTEGER :: n, n1x, t
    REAL*8 :: grad
    
    n1x = SIZE(x1)
  
    !
    ! ... locate the grid points near the topographic points
    ! ... and interpolate linearly the profile
    !
    t = 1

    DO n = 1, n1x

       IF (x1(n) <= x2) t = n

    END DO
    
    IF (t==1 .OR. t==n1x) THEN

       f2 = f1(t)

    ELSE

       grad = (f1(t+1)-f1(t))/(x1(t+1)-x1(t))
       f2 = f1(t) + (x2-x1(t)) * grad

    END IF

    RETURN
    
  END SUBROUTINE interp_1d_scalar
  
  !------------------------------------------------------------------

END MODULE meteo_module
!----------------------------------------------------------------------