meteo.f90

Go to the documentation of this file.

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

MODULE meteo_module

  USE variables, ONLY: gi   ! Grav acceleration 
  USE variables, ONLY : wp

  IMPLICIT NONE
    
  REAL(wp) :: h0    
  REAL(wp) :: h1    
  REAL(wp) :: h2    
  REAL(wp) :: h3    
  REAL(wp) :: h4    
  REAL(wp) :: h5    
  REAL(wp) :: h6    

  REAL(wp) :: u0    
  REAL(wp) :: u1    
  REAL(wp) :: u2    
  REAL(wp) :: u3    
  REAL(wp) :: u4    
  REAL(wp) :: u5    
  REAL(wp) :: u6    
  
  REAL(wp) :: sp_hu0, sp_hu1, sp_hu2, sp_hu3, sp_hu4, sp_hu5, sp_hu6 

  REAL(wp) :: rel_hu
  
  REAL(wp) :: gamma_m

  REAL(wp) :: gamma_d

  REAL(wp) :: rho_dry
  
  REAL(wp) :: rh
  
  REAL(wp) :: cos_theta , sin_theta

  REAL(wp) :: rho_atm0

  REAL(wp) :: u_atm   

  REAL(wp) :: u_wind
  REAL(wp) :: v_wind
  
  REAL(wp) :: rho_atm  

  REAL(wp) :: sphu_atm

  REAL(wp) :: sphu_atm0

  
  REAL(wp) :: visc_atm 

  REAL(wp) :: visc_atm0

  REAL(wp) :: ta      

  REAL(wp) :: pa      

  REAL(wp) :: dpdz     

  REAL(wp) :: dtdz     

  REAL(wp) :: rair

  REAL(wp) :: cpair
  
  REAL(wp), PARAMETER :: t_ref = 273.15_wp
  
  REAL(wp), PARAMETER :: h_wv0 = 2.501d6
  
  REAL(wp), PARAMETER :: c_wv = 1996.0_wp
  
  REAL(wp), PARAMETER :: h_lw0 = 3.337d5
  
  REAL(wp), PARAMETER :: c_lw = 4187.0_wp

  REAL(wp), PARAMETER :: c_ice = 2108.0_wp
  
  REAL(wp), PARAMETER :: da_mol_wt = 0.029_wp
  
  REAL(wp), PARAMETER :: wv_mol_wt = 0.018_wp
   
  REAL(wp), PARAMETER :: rwv = 462

  INTEGER :: n_atm_profile

  REAL(wp), ALLOCATABLE :: atm_profile(:,:)

  CHARACTER*10 :: read_atm_profile

  REAL(wp) :: u_max , z_r , exp_wind

  REAL(wp), ALLOCATABLE :: rho_atm_month_lat(:) , pres_atm_month_lat(:) ,   &
       temp_atm_month_lat(:) , temp_atm_month(:,:)

  REAL(wp), ALLOCATABLE :: h_levels(:)
  
  REAL(wp) :: wind_mult_coeff

  SAVE

CONTAINS

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

  SUBROUTINE initialize_meteo

    IMPLICIT NONE


    IF ( read_atm_profile .EQ. 'standard' ) THEN

       ! tropospere base altitude (m)
       h0 = 0.0_wp 
       u0 = 0.0_wp
       sp_hu0 = sphu_atm0
       
       ! ------------------------ atmosphere layer 1 ----------------------------

       ! tropopause base altitude
       h1 = 11000.0_wp
       ! specific humidity (kg/kg) at top of layer 1
       sp_hu1 = 2.0e-6_wp
       u1 = u_max
       
       ! ------------------------ atmosphere layer 2 ----------------------------
       h2 = 20000.0_wp
       ! specific humidity (kg/kg) at top of layer 2
       sp_hu2 = 2.6e-6_wp
       u2 = 10.0_wp

       ! ------------------------ atmosphere layer 3 ----------------------------
       h3 = 32000.0_wp
       ! specific humidity (kg/kg) at top of layer 3
       sp_hu3 = 3.2e-6_wp

       ! ------------------------ atmosphere layer 4 ----------------------------
       h4 = 47000.0_wp
       ! specific humidity (kg/kg) at top of layer 4
       sp_hu4 = 3.2e-6_wp
        
       ! ------------------------ atmosphere layer 5 ----------------------------
       h5 = 51000.0_wp
       ! specific humidity (kg/kg) at top of layer 5
       sp_hu5 = 3.2e-6_wp

       ! ------------------------ atmosphere layer 6 ----------------------------
       h6 = 71000.0_wp
       ! specific humidity (kg/kg) at top of layer 6
       sp_hu6 = 2.4e-6_wp
       ! Ref: WINDS AT ALTITUDES UP TO 80 KILOMETERS (Fig. 2)
       u6 = 65.0_wp

    END IF

    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(wp) :: const, const1, t1, p1, p2
    REAL(wp) :: const2

    REAL(wp) :: WE_wind , NS_wind

    REAL(wp) :: Cs

    ! Density correction factor
    REAL(wp) :: K

    REAL(wp) :: sp_hu_avg
    
    REAL(wp) :: hu_ratio

    REAL(wp) :: h_bot , h_top
    REAL(wp) :: sphu_bot , sphu_top

    REAL(wp) :: T_ref , t0
    REAL(wp) :: el , es , e_sl

    REAL(wp) :: p_wv , p_da
    REAL(wp) :: n_wv , n_da
    REAL(wp) :: x_wv , x_da

    
    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 specific humifity profile
       CALL interp_1d_scalar(atm_profile(1,:), atm_profile(5,:), z, sphu_atm)

       ! 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_wind = we_wind
       v_wind = -ns_wind
       
       IF ( ( we_wind .EQ. 0.0_wp ) .AND. ( ns_wind .EQ. 0.0_wp ) ) THEN

          we_wind = 1.e-15_wp
          ns_wind = 1.e-15_wp
          
       END IF
       
       u_atm = sqrt( 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 - vent_height ) .LE. z_r ) THEN

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

       ELSE

          u_atm = u_max

       END IF
       
       cos_theta = 1.0_wp
       sin_theta = 0.0_wp

    ELSEIF ( read_atm_profile .EQ. 'standard' ) THEN

       IF ( z <= h1 ) THEN

          ! ... Troposphere
          u_atm = u0 + (u1-u0) * (z-h0) / (h1-h0)

          h_bot = h0
          h_top = h1
          sphu_bot = sp_hu0
          sphu_top = sp_hu1
          
          gamma_d = 6.5e-3_wp

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

          ! ... Tropopause
          u_atm = u1 + (u2-u1) * (z-h1) / (h2-h1)

          h_bot = h1
          h_top = h2
          sphu_bot = sp_hu1
          sphu_top = sp_hu2

          gamma_d = 0.0_wp

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

          ! ... Stratosphere
          u_atm = u2 + (u6-u2) * (z-h2) / (h6-h2)

          h_bot = h2
          h_top = h3
          sphu_bot = sp_hu2
          sphu_top = sp_hu3

          gamma_d = -1.0e-3_wp

       ELSE IF (z > h3 .AND. z <= h4) THEN

          u_atm = u2 + (u6-u2) * (z-h2) / (h6-h2);

          h_bot = h3
          h_top = h4
          sphu_bot = sp_hu3
          sphu_top = sp_hu4
          
          gamma_d = -2.8e-3_wp
          
       ELSE IF (z > h4 .AND. z <= h5) THEN

          u_atm = u2 + (u6-u2) * (z-h2) / (h6-h2)

          h_bot = h4
          h_top = h5
          sphu_bot = sp_hu4
          sphu_top = sp_hu5

          gamma_d = 0.0_wp
          
       ELSE IF (z > h5 .AND. z <= h6) THEN

          u_atm = u2 + (u6-u2) * (z-h2) / (h6-h2);

          h_bot = h5
          h_top = h6
          sphu_bot = sp_hu5
          sphu_top = sp_hu6

          gamma_d = +2.8e-3_wp

       ENDIF

       IF ( sphu_atm0 .GT. 0.0_wp ) THEN
       
          ! log of specific humidity is assumed to vary linearly
          sphu_atm = exp( log(sphu_bot) + ( log(sphu_top) - log(sphu_bot) ) *   &
               (z-h_bot) / (h_top-h_bot) )

       ELSE

          t_ref = 273.16_wp
          t0 = t_ref - 40.0_wp

          ! saturation pressure of vapor over liquid Pa
          el = 611.2_wp * exp( 17.67_wp * ( ta-273.16_wp ) / ( ta - 29.65_wp ) )
          
          ! saturation pressure of vapor over ice Pa
          es = -9.097_wp * ( (273.16_wp / ta ) - 1.0_wp ) - 3.566_wp *          &
               log10(273.16_wp / ta) + 0.876_wp * ( 1.0_wp - (ta / 273.16_wp) )
          
          es = 611.22_wp * ( 10.0_wp**es ) 
                    
          ! saturation pressure of vapor
          IF ( ta .GE. t_ref ) THEN
             
             e_sl = el
             
          ELSEIF ( ta .LE. t0 ) THEN
             
             e_sl = es
             
          ELSE
             
             e_sl = es + ( ta - t0 ) / ( t_ref - t0 ) * ( el - es )
             
          END IF
          
          p_wv = min(pa, rel_hu*e_sl)
          ! dry air partial pressure
          p_da = pa - p_wv  
          !molar fraction of water vapor
          n_wv = p_wv / pa
          !molar fraction dry air
          n_da = p_da / pa
          !mass fraction of da 
          x_da = (n_da * da_mol_wt) / (n_da * da_mol_wt + n_wv * wv_mol_wt )
          !mass fraction of wv
          x_wv = (n_wv * wv_mol_wt) / (n_da * da_mol_wt + n_wv * wv_mol_wt )
          !specific humidity
          sphu_atm = x_wv / (x_da+x_wv)

       END IF

       ! WRITE(*,*) sphu_atm 
       
       ! Density of dry air
       rho_dry = pa / ( rair*ta )

       ! Density correction factor
       k = 1.0_wp / ( 1.0_wp + ( rwv/rair - 1.0_wp ) * sphu_atm )

       ! Density of mixure (dry air and water vapor)
       rho_atm = k * rho_dry

       ! Lapse rate corrected for specific humidity
       gamma_m = gamma_d * (1.0_wp - 0.856_wp * sphu_atm )
       
       u_wind = u_atm  
       v_wind = 0.0_wp
       
       IF ( u_max .GT. 0.0_wp ) THEN
          
          cos_theta = u_wind/u_max 
          sin_theta = v_wind/u_max
          
       ELSE
          
          cos_theta = 1.0_wp
          sin_theta = 0.0_wp
          
       END IF
              
    END IF

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

    IF ( verbose_level .GE. 1 ) THEN
       
       WRITE(*,*) 'Height (asl) = ',z
       WRITE(*,*) 'Ambient temperature (K) = ',ta
       WRITE(*,*) 'Ambient pressure (Pa) = ', pa
       WRITE(*,*) 'Dry atmosphere density (kg m-3) = ',rho_dry
       WRITE(*,*) 'Moist atmosphere density (kg m-3) = ',rho_atm
       WRITE(*,*) 'Atmosphere viscosity = ',visc_atm
       WRITE(*,*) 'Wind speed (m s-1) = ',u_atm
       READ(*,*)

    END IF

    RETURN

  END SUBROUTINE zmet

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

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

    DO n = 1, n1x

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

    END DO
    
    IF ( t.EQ.0 ) THEN

       f2 = f1(1)

    ELSEIF ( t.EQ.n1x ) THEN

       f2 = f1(n1x)

    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
!----------------------------------------------------------------------