Ignore:
Timestamp:
Jan 27, 2016, 10:42:32 AM (8 years ago)
Author:
idelkadi
Message:

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

  • mise a jour des sources pour ISCCP, CALIPSO et PARASOL
  • prise en compte des changements de phases pour les nuages (Calipso)
  • rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ5/trunk/libf/phylmd/cosp/radar_simulator.F90

    r1907 r2428  
    1   subroutine radar_simulator(freq,k2,do_ray,use_gas_abs,use_mie_table,mt, &
    2     nhclass,hp,nprof,ngate,nsizes,D,hgt_matrix,hm_matrix,re_matrix,p_matrix,t_matrix, &
    3     rh_matrix,Ze_non,Ze_ray,h_atten_to_vol,g_atten_to_vol,dBZe, &
     1  subroutine radar_simulator( &
     2    hp, &
     3    nprof,ngate, &
     4    undef, &
     5    hgt_matrix,hm_matrix,re_matrix,Np_matrix, &
     6    p_matrix,t_matrix,rh_matrix, &
     7    Ze_non,Ze_ray,a_to_vol,g_to_vol,dBZe, &
    48    g_to_vol_in,g_to_vol_out)
    5 
    6 !     rh_matrix,Ze_non,Ze_ray,kr_matrix,g_atten_to_vol,dBZe)
    79 
    810  use m_mrgrnk
    … …  
    1113  use optics_lib
    1214  use radar_simulator_types
     15  use scale_LUTs_io
    1316  implicit none
    1417 
    1518! Purpose:
     19!
    1620!   Simulates a vertical profile of radar reflectivity
    17 !   Part of QuickBeam v1.04 by John Haynes & Roger Marchand
     21!   Originally Part of QuickBeam v1.04 by John Haynes & Roger Marchand.
     22!   but has been substantially modified since that time by
     23!   Laura Fowler and Roger Marchand (see modifications below).
    1824!
    1925! Inputs:
    20 !   [freq]            radar frequency (GHz), can be anything unless
    21 !                     use_mie_table=1, in which case one of 94,35,13.8,9.6,3
    22 !   [k2]              |K|^2, the dielectric constant, set to -1 to use the
    23 !                     frequency dependent default
    24 !   [do_ray]          1=do Rayleigh calcs, 0=not
    25 !   [use_gas_abs]     1=do gaseous abs calcs, 0=not,
    26 !                     2=use same as first profile (undocumented)
    27 !   [use_mie_table]   1=use Mie tables, 0=not
    28 !   [mt]              Mie look up table
    29 !   [nhclass]         number of hydrometeor types
    30 !   [hp]              structure that defines hydrometeor types
     26!
     27!   [hp]              structure that defines hydrometeor types and other radar properties
     28!
    3129!   [nprof]           number of hydrometeor profiles
    3230!   [ngate]           number of vertical layers
    33 !   [nsizes]          number of discrete particles in [D]
    34 !   [D]               array of discrete particles (um)
    35 !
     31!
     32!   [undef]           missing data value
    3633!   (The following 5 arrays must be in order from closest to the radar
    3734!    to farthest...)
     35!
    3836!   [hgt_matrix]      height of hydrometeors (km)
     37!   [p_matrix]        pressure profile (hPa)
     38!   [t_matrix]        temperature profile (K)
     39!   [rh_matrix]       relative humidity profile (%) -- only needed if gaseous aborption calculated.
     40!
    3941!   [hm_matrix]       table of hydrometeor mixing rations (g/kg)
    40 !   [re_matrix]       OPTIONAL table of hydrometeor effective radii (microns)
    41 !   [p_matrix]        pressure profile (hPa)
    42 !   [t_matrix]        temperature profile (C)
    43 !   [rh_matrix]       relative humidity profile (%)
     42!   [re_matrix]       table of hydrometeor effective radii.  0 ==> use defaults. (units=microns)   
     43!   [Np_matrix]       table of hydrometeor number concentration.  0 ==> use defaults. (units = 1/kg)
    4444!
    4545! Outputs:
     46!
    4647!   [Ze_non]          radar reflectivity without attenuation (dBZ)
    4748!   [Ze_ray]          Rayleigh reflectivity (dBZ)
    … …  
    5960! Created:
    6061!   11/28/2005  John Haynes (haynes@atmos.colostate.edu)
     62!
    6163! Modified:
    62 !   09/2006  placed into subroutine form, scaling factors (Roger Marchand,JMH)
     64!   09/2006  placed into subroutine form (Roger Marchand,JMH)
    6365!   08/2007  added equivalent volume spheres, Z and N scalling most distrubtion types (Roger Marchand)
    6466!   01/2008  'Do while' to determine if hydrometeor(s) present in volume
    6567!             changed for vectorization purposes (A. Bodas-Salcedo)
    66 
     68!
     69!   07/2010  V3.0 ... Modified to load or save scale factors to disk as a Look-Up Table (LUT)
     70!  ... All hydrometeor and radar simulator properties now included in hp structure
     71!  ... hp structure should be initialized by call to radar_simulator_init prior
     72!  ... to calling this subroutine. 
     73!     Also ... Support of Morrison 2-moment style microphyscis (Np_matrix) added
     74!  ... Changes implement by Roj Marchand following work by Laura Fowler
     75!
     76!   10/2011  Modified ngate loop to go in either direction depending on flag
     77!     hp%radar_at_layer_one.  This affects the direction in which attenuation is summed.
     78!
     79!     Also removed called to AVINT for gas and hydrometeor attenuation and replaced with simple
     80!     summation. (Roger Marchand)
     81!
     82!
    6783! ----- INPUTS ----- 
    68   type(mie), intent(in) :: mt
     84 
     85  logical, parameter  ::  DO_LUT_TEST = .false.
     86  logical, parameter  ::  DO_NP_TEST = .false.
     87
    6988  type(class_param), intent(inout) :: hp
    70   real*8, intent(in) :: freq,k2
    71   integer, intent(in) ::  do_ray,use_gas_abs,use_mie_table, &
    72     nhclass,nprof,ngate,nsizes
    73   real*8, dimension(nsizes), intent(in) :: D
    74   real*8, dimension(nprof,ngate), intent(in) :: hgt_matrix, p_matrix, &
    75     t_matrix,rh_matrix
    76   real*8, dimension(nhclass,nprof,ngate), intent(in) :: hm_matrix
    77   real*8, dimension(nhclass,nprof,ngate), intent(inout) :: re_matrix
    78    
     89
     90  integer, intent(in) ::  nprof,ngate
     91
     92  real undef
     93  real*8, dimension(nprof,ngate), intent(in) :: &
     94    hgt_matrix, p_matrix,t_matrix,rh_matrix
     95
     96  real*8, dimension(hp%nhclass,nprof,ngate), intent(in) :: hm_matrix
     97  real*8, dimension(hp%nhclass,nprof,ngate), intent(inout) :: re_matrix
     98  real*8, dimension(hp%nhclass,nprof,ngate), intent(in)    :: Np_matrix
     99
    79100! ----- OUTPUTS -----
    80101  real*8, dimension(nprof,ngate), intent(out) :: Ze_non,Ze_ray, &
    81         g_atten_to_vol,dBZe,h_atten_to_vol
     102       g_to_vol,dBZe,a_to_vol
    82103
    83104! ----- OPTIONAL -----
    84   real*8, optional, dimension(ngate,nprof) :: &
     105  real*8, optional, dimension(nprof,ngate) :: &
    85106  g_to_vol_in,g_to_vol_out ! integrated atten due to gases, r>v (dB). This allows to output and then input
    86107                           ! the same gaseous absorption in different calls. Optional to allow compatibility
    87108                           ! with original version. A. Bodas April 2008.
    88        
    89109!  real*8, dimension(nprof,ngate) :: kr_matrix
    90110
    91111! ----- INTERNAL -----
     112
     113  real, parameter :: one_third = 1.0/3.0
     114  real*8 :: t_kelvin
    92115  integer :: &
    93   phase, &                      ! 0=liquid, 1=ice
    94   ns                            ! number of discrete drop sizes
    95 
    96   integer*4, dimension(ngate) :: &
    97   hydro                         ! 1=hydrometeor in vol, 0=none
     116  phase, & ! 0=liquid, 1=ice
     117  ns       ! number of discrete drop sizes
     118
     119  logical :: hydro      ! true=hydrometeor in vol, false=none
    98120  real*8 :: &
    99   rho_a, &                      ! air density (kg m^-3)
    100   gases                         ! function: 2-way gas atten (dB/km)
     121  rho_a, &   ! air density (kg m^-3)
     122  gases      ! function: 2-way gas atten (dB/km)
    101123
    102124  real*8, dimension(:), allocatable :: &
    103   Di, Deq, &                    ! discrete drop sizes (um)
    104   Ni, Ntemp, &                  ! discrete concentrations (cm^-3 um^-1)
    105   rhoi                          ! discrete densities (kg m^-3)
    106  
    107   real*8, dimension(ngate) :: &
    108   z_vol, &                      ! effective reflectivity factor (mm^6/m^3)
     125  Di, Deq, &   ! discrete drop sizes (um)
     126  Ni, &        ! discrete concentrations (cm^-3 um^-1)
     127  rhoi         ! discrete densities (kg m^-3)
     128
     129  real*8, dimension(nprof, ngate) :: &
     130  z_vol, &      ! effective reflectivity factor (mm^6/m^3)
    109131  z_ray, &                      ! reflectivity factor, Rayleigh only (mm^6/m^3)
    110   kr_vol, &                     ! attenuation coefficient hydro (dB/km)
    111   g_vol, &                      ! attenuation coefficient gases (dB/km)
    112   a_to_vol, &                   ! integrated atten due to hydometeors, r>v (dB)
    113   g_to_vol                      ! integrated atten due to gases, r>v (dB)
    114    
    115  
     132  kr_vol, &     ! attenuation coefficient hydro (dB/km)
     133  g_vol         ! attenuation coefficient gases (dB/km)
     134
     135
    116136  integer,parameter :: KR8 = selected_real_kind(15,300)
    117137  real*8, parameter :: xx = -1.0_KR8
    118138  real*8,  dimension(:), allocatable :: xxa
    119   real*8 :: kr, ze, zr, pi, scale_factor, tc, Re, ld, tmp1, ze2, kr2,apm,bpm
     139  real*8 :: kr, ze, zr, pi, scale_factor, tc, Re, ld, tmp1, ze2, kr2, apm, bpm
     140  real*8 :: half_a_atten_current,half_a_atten_above
     141  real*8 :: half_g_atten_current,half_g_atten_above
    120142  integer*4 :: tp, i, j, k, pr, itt, iff
    121143
    122   real*8 bin_length,step,base,step_list(25),base_list(25)
     144  real*8    step,base, Np
    123145  integer*4 iRe_type,n,max_bin
    124  
     146
     147  integer   start_gate,end_gate,d_gate
     148
    125149  logical :: g_to_vol_in_present, g_to_vol_out_present
    126        
     150
    127151  ! Logicals to avoid calling present within the loops
    128152  g_to_vol_in_present  = present(g_to_vol_in)
    129153  g_to_vol_out_present = present(g_to_vol_out)
    130  
    131     ! set up Re bins for z_scalling
    132         bin_length=50;
    133         max_bin=25
    134 
    135         step_list(1)=1
    136         base_list(1)=75
    137         do j=2,max_bin
    138                 step_list(j)=3*(j-1);
    139                 if(step_list(j)>bin_length) then
    140                         step_list(j)=bin_length;
    141                 endif
    142                 base_list(j)=base_list(j-1)+floor(bin_length/step_list(j-1));
    143         enddo
    144 
     154
     155  !
     156  ! load scaling matricies from disk -- but only the first time this subroutine is called
     157  !
     158  if(hp%load_scale_LUTs) then
     159    call load_scale_LUTs(hp)
     160    hp%load_scale_LUTs=.false.
     161    hp%Z_scale_added_flag = .false. ! will be set true if scaling Look Up Tables are modified during run
     162  endif
    145163
    146164  pi = acos(-1.0)
    147   if (use_mie_table == 1) iff = infind(mt%freq,freq,sort=1)
    148 
    149        
     165
     166!   ----- Initialisation -----
     167  g_to_vol = 0.0
     168  a_to_vol = 0.0
     169  z_vol    = 0.0
     170  z_ray    = 0.0
     171  kr_vol   = 0.0
     172
     173!   // loop over each range gate (ngate) ... starting with layer closest to the radar !
     174  if(hp%radar_at_layer_one) then
     175    start_gate=1
     176    end_gate=ngate
     177    d_gate=1
     178  else
     179    start_gate=ngate
     180    end_gate=1
     181    d_gate=-1
     182  endif
     183  do k=start_gate,end_gate,d_gate
    150184  ! // loop over each profile (nprof)
    151   do pr=1,nprof
    152 
    153 !   ----- calculations for each volume -----
    154     z_vol(:) = 0
    155     z_ray(:) = 0
    156     kr_vol(:) = 0
    157     hydro(:) = 0   
    158 
    159 !   // loop over eacho range gate (ngate)
    160     do k=1,ngate
    161  
     185    do pr=1,nprof
     186      t_kelvin = t_matrix(pr,k)
    162187!     :: determine if hydrometeor(s) present in volume
    163       hydro(k) = 0
    164       do j=1,nhclass ! Do while changed for vectorization purposes (A. B-S)
     188      hydro = .false.
     189      do j=1,hp%nhclass
    165190        if ((hm_matrix(j,pr,k) > 1E-12) .and. (hp%dtype(j) > 0)) then
    166           hydro(k) = 1
     191          hydro = .true.
    167192          exit
    168193        endif
    169194      enddo
    170195
    171       if (hydro(k) == 1) then
    172 !     :: if there is hydrometeor in the volume           
    173 
    174         rho_a = (p_matrix(pr,k)*100.)/(287*(t_matrix(pr,k)+273.15))
    175 
     196!     :: if there is hydrometeor in the volume
     197      if (hydro) then
     198
     199        rho_a = (p_matrix(pr,k)*100.)/(287.0*(t_kelvin))
    176200!       :: loop over hydrometeor type
    177         do tp=1,nhclass
    178 
     201        do tp=1,hp%nhclass
    179202          if (hm_matrix(tp,pr,k) <= 1E-12) cycle
    180 
    181           phase = hp%phase(tp)
    182           if(phase==0) then
    183                 itt = infind(mt_ttl,t_matrix(pr,k))
    184           else
    185                 itt = infind(mt_tti,t_matrix(pr,k))
     203          phase = hp%phase(tp)
     204          if (phase==0) then
     205            itt = infind(hp%mt_ttl,t_kelvin)
     206          else
     207            itt = infind(hp%mt_tti,t_kelvin)
     208          endif
     209          if (re_matrix(tp,pr,k).eq.0) then
     210            call calc_Re(hm_matrix(tp,pr,k),Np_matrix(tp,pr,k),rho_a, &
     211              hp%dtype(tp),hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), &
     212              hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp),Re)
     213            re_matrix(tp,pr,k)=Re
     214          else
     215            if (Np_matrix(tp,pr,k)>0) then
     216               print *, 'Warning: Re and Np set for the same ', &
     217                        'volume & hydrometeor type.  Np is being ignored.'
     218            endif
     219            Re = re_matrix(tp,pr,k)
     220          endif
     221
     222          iRe_type=1
     223          if(Re.gt.0) then
     224            ! determine index in to scale LUT
     225            !
     226            ! distance between Re points (defined by "base" and "step") for
     227            ! each interval of size Re_BIN_LENGTH
     228            ! Integer asignment, avoids calling floor intrinsic
     229            n=Re/Re_BIN_LENGTH
     230            if (n>=Re_MAX_BIN) n=Re_MAX_BIN-1
     231            step=hp%step_list(n+1)
     232            base=hp%base_list(n+1)
     233            iRe_type=Re/step
     234            if (iRe_type.lt.1) iRe_type=1
     235
     236            Re=step*(iRe_type+0.5)      ! set value of Re to closest value allowed in LUT.
     237            iRe_type=iRe_type+base-int(n*Re_BIN_LENGTH/step)
     238
     239            ! make sure iRe_type is within bounds
     240            if (iRe_type.ge.nRe_types) then
     241!               write(*,*) 'Warning: size of Re exceed value permitted ', &
     242!                    'in Look-Up Table (LUT).  Will calculate. '
     243               ! no scaling allowed
     244               iRe_type=nRe_types
     245               hp%Z_scale_flag(tp,itt,iRe_type)=.false.
     246            else
     247               ! set value in re_matrix to closest values in LUT
     248              if (.not. DO_LUT_TEST) re_matrix(tp,pr,k)=Re
     249            endif
     250          endif
     251          ! use Ze_scaled, Zr_scaled, and kr_scaled ... if know them
     252          ! if not we will calculate Ze, Zr, and Kr from the distribution parameters
     253          if( (.not. hp%Z_scale_flag(tp,itt,iRe_type)) .or. DO_LUT_TEST)  then
     254            ! :: create a distribution of hydrometeors within volume
     255            select case(hp%dtype(tp))
     256              case(4)
     257                ns = 1
     258                allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
     259                Di = hp%p1(tp)
     260                Ni = 0.
     261              case default
     262                ns = nd   ! constant defined in radar_simulator_types.f90
     263                allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
     264                Di = hp%D
     265                Ni = 0.
     266            end select
     267            call dsd(hm_matrix(tp,pr,k),re_matrix(tp,pr,k),Np_matrix(tp,pr,k), &
     268                     Di,Ni,ns,hp%dtype(tp),rho_a,t_kelvin, &
     269                     hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), &
     270                     hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp))
     271
     272            ! calculate particle density
     273            if (phase == 1) then
     274              if (hp%rho(tp) < 0) then
     275                ! Use equivalent volume spheres.
     276                hp%rho_eff(tp,1:ns,iRe_type) = 917                              ! solid ice == equivalent volume approach
     277                Deq = ( ( 6/pi*hp%apm(tp)/917 ) ** (1.0/3.0) ) * ( (Di*1E-6) ** (hp%bpm(tp)/3.0) )  * 1E6
     278                ! alternative is to comment out above two lines and use the following block
     279                ! MG Mie approach - adjust density of sphere with D = D_characteristic to match particle density
     280                !
     281                ! hp%rho_eff(tp,1:ns,iRe_type) = (6/pi)*hp%apm(tp)*(Di*1E-6)**(hp%bpm(tp)-3)   !MG Mie approach
     282
     283                ! as the particle size gets small it is possible that the mass to size relationship of
     284                ! (given by power law in hclass.data) can produce impossible results
     285                ! where the mass is larger than a solid sphere of ice. 
     286                ! This loop ensures that no ice particle can have more mass/density larger than an ice sphere.
     287                ! do i=1,ns
     288                ! if(hp%rho_eff(tp,i,iRe_type) > 917 ) then
     289                ! hp%rho_eff(tp,i,iRe_type) = 917
     290                ! endif
     291                ! enddo
     292              else
     293                ! Equivalent volume sphere (solid ice rho_ice=917 kg/m^3).
     294                hp%rho_eff(tp,1:ns,iRe_type) = 917
     295                Deq=Di * ((hp%rho(tp)/917)**(1.0/3.0))
     296                ! alternative ... coment out above two lines and use the following for MG-Mie
     297                ! hp%rho_eff(tp,1:ns,iRe_type) = hp%rho(tp)   !MG Mie approach
     298              endif
     299            else
     300              ! I assume here that water phase droplets are spheres.
     301              ! hp%rho should be ~ 1000  or hp%apm=524 .and. hp%bpm=3
     302              Deq = Di
     303            endif
     304
     305            ! calculate effective reflectivity factor of volume
     306            xxa = -9.9
     307            rhoi = hp%rho_eff(tp,1:ns,iRe_type)
     308            call zeff(hp%freq,Deq,Ni,ns,hp%k2,t_kelvin,phase,hp%do_ray, &
     309                      ze,zr,kr,xxa,xxa,rhoi)
     310
     311            ! test code ... compare Np value input to routine with sum of DSD
     312            ! NOTE: if .not. DO_LUT_TEST, then you are checking the LUT approximation
     313            ! not just the DSD representation given by Ni
     314            if(Np_matrix(tp,pr,k)>0 .and. DO_NP_TEST ) then
     315              Np = path_integral(Ni,Di,1,ns-1)/rho_a*1E6
     316              ! Note: Representation is not great or small Re < 2
     317              if( (Np_matrix(tp,pr,k)-Np)/Np_matrix(tp,pr,k)>0.1 ) then
     318                write(*,*) 'Error: Np input does not match sum(N)'
     319                write(*,*) tp,pr,k,Re,Ni(1),Ni(ns),10*log10(ze)
     320                write(*,*) Np_matrix(tp,pr,k),Np,(Np_matrix(tp,pr,k)-Np)/Np_matrix(tp,pr,k)
     321                write(*,*)
     322              endif
     323            endif
     324
     325            deallocate(Di,Ni,rhoi,xxa,Deq)
     326
     327            ! LUT test code
     328            ! This segment of code compares full calculation to scaling result
     329            if ( hp%Z_scale_flag(tp,itt,iRe_type) .and. DO_LUT_TEST )  then
     330              scale_factor=rho_a*hm_matrix(tp,pr,k)
     331              ! if more than 2 dBZe difference print error message/parameters.
     332              if ( abs(10*log10(ze) - 10*log10(hp%Ze_scaled(tp,itt,iRe_type) * &
     333                   scale_factor)) > 2 ) then
     334                write(*,*) 'Roj Error: ',tp,itt,iRe_type,hp%Z_scale_flag(tp,itt,iRe_type),n,step,base
     335                write(*,*) 10*log10(ze),10*log10(hp%Ze_scaled(tp,itt,iRe_type) * scale_factor)
     336                write(*,*) hp%Ze_scaled(tp,itt,iRe_type),scale_factor
     337                write(*,*) re_matrix(tp,pr,k),Re
     338                write(*,*)
     339              endif
     340            endif
     341
     342          else ! can use z scaling
     343            scale_factor=rho_a*hm_matrix(tp,pr,k)
     344            zr = hp%Zr_scaled(tp,itt,iRe_type) * scale_factor
     345            ze = hp%Ze_scaled(tp,itt,iRe_type) * scale_factor
     346            kr = hp%kr_scaled(tp,itt,iRe_type) * scale_factor
     347          endif  ! end z_scaling
     348
     349          kr_vol(pr,k) = kr_vol(pr,k) + kr
     350          z_vol(pr,k)  = z_vol(pr,k)  + ze
     351          z_ray(pr,k)  = z_ray(pr,k)  + zr
     352
     353          ! construct Ze_scaled, Zr_scaled, and kr_scaled ... if we can
     354          if ( .not. hp%Z_scale_flag(tp,itt,iRe_type) ) then
     355            if (iRe_type>1) then
     356              scale_factor=rho_a*hm_matrix(tp,pr,k)
     357              hp%Ze_scaled(tp,itt,iRe_type) = ze/ scale_factor
     358              hp%Zr_scaled(tp,itt,iRe_type) = zr/ scale_factor
     359              hp%kr_scaled(tp,itt,iRe_type) = kr/ scale_factor
     360              hp%Z_scale_flag(tp,itt,iRe_type) = .true.
     361              hp%Z_scale_added_flag(tp,itt,iRe_type)=.true.
     362            endif
     363          endif
     364
     365        enddo   ! end loop of tp (hydrometeor type)
     366
     367      else
     368!     :: volume is hydrometeor-free     
     369        kr_vol(pr,k) = 0
     370        z_vol(pr,k)  = undef
     371        z_ray(pr,k)  = undef
    186372      endif
    187373
    188           ! calculate Re if we have an exponential distribution with fixed No ... precipitation type particle
    189           if( hp%dtype(tp)==2 .and. abs(hp%p2(tp)+1) < 1E-8)  then
    190 
    191                 apm=hp%apm(tp)
    192                 bpm=hp%bpm(tp)
    193 
    194                 if ((hp%rho(tp) > 0) .and. (apm < 0)) then
    195                         apm = (pi/6)*hp%rho(tp)
    196                         bpm = 3.
    197                 endif
    198 
    199                 tmp1 = 1./(1.+bpm)
    200                 ld = ((apm*gamma(1.+bpm)*hp%p1(tp))/(rho_a*hm_matrix(tp,pr,k)*1E-3))**tmp1
    201                
    202                 Re = 1.5E6/ld
    203                
    204                 re_matrix(tp,pr,k) = Re;
    205 
    206           endif
    207  
    208           if(re_matrix(tp,pr,k).eq.0) then
    209 
    210                 iRe_type=1
    211                 Re=0
    212           else
    213                 iRe_type=1
    214                 Re=re_matrix(tp,pr,k)
    215                
    216                 n=floor(Re/bin_length)
    217                 if(n==0) then
    218                         if(Re<25) then
    219                                 step=0.5
    220                                 base=0
    221                         else                   
    222                                 step=1
    223                                 base=25
    224                         endif
    225                 else
    226                         if(n>max_bin) then
    227                                 n=max_bin       
    228                         endif
    229 
    230                         step=step_list(n)
    231                         base=base_list(n)
    232                 endif
    233 
    234                 iRe_type=floor(Re/step)
    235 
    236                 if(iRe_type.lt.1) then 
    237                         iRe_type=1                     
    238                 endif
    239 
    240                 Re=step*(iRe_type+0.5)
    241                 iRe_type=iRe_type+base-floor(n*bin_length/step)
    242 
    243                 ! make sure iRe_type is within bounds
    244                 if(iRe_type.ge.nRe_types) then 
    245 
    246                         ! print *, tp, re_matrix(tp,pr,k), Re, iRe_type
    247 
    248                         ! no scaling allowed
    249                         Re=re_matrix(tp,pr,k)
    250 
    251                         iRe_type=nRe_types
    252                         hp%z_flag(tp,itt,iRe_type)=.false.
    253                         hp%scaled(tp,iRe_type)=.false.                 
    254                 endif
    255           endif
    256        
    257           ! use Ze_scaled, Zr_scaled, and kr_scaled ... if know them
    258           ! if not we will calculate Ze, Zr, and Kr from the distribution parameters
    259           if( .not. hp%z_flag(tp,itt,iRe_type) )  then
    260          
    261 !         :: create a distribution of hydrometeors within volume         
    262           select case(hp%dtype(tp))
    263           case(4)
    264             ns = 1
    265             allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
    266             if (use_mie_table == 1) allocate(mt_qext(ns),mt_qbsca(ns),Ntemp(ns))
    267             Di = hp%p1(tp)
    268             Ni = 0.
    269           case default
    270             ns = nsizes           
    271             allocate(Di(ns),Ni(ns),rhoi(ns),xxa(ns),Deq(ns))
    272             if (use_mie_table == 1) allocate(mt_qext(ns),mt_qbsca(ns),Ntemp(ns))           
    273             Di = D
    274             Ni = 0.
    275           end select
    276 
    277 !         :: create a DSD (using scaling factor if applicable)
    278           ! hp%scaled(tp,iRe_type)=.false.   ! turn off N scaling
    279 
    280           call dsd(hm_matrix(tp,pr,k),Re,Di,Ni,ns,hp%dtype(tp),rho_a, &
    281             t_matrix(pr,k),hp%dmin(tp),hp%dmax(tp),hp%apm(tp),hp%bpm(tp), &
    282             hp%rho(tp),hp%p1(tp),hp%p2(tp),hp%p3(tp),hp%fc(tp,1:ns,iRe_type), &
    283             hp%scaled(tp,iRe_type))
    284 
    285 !         :: calculate particle density
    286           ! if ((hp%rho_eff(tp,1,iRe_type) < 0) .and. (phase == 1)) then
    287           if (phase == 1) then
    288             if (hp%rho(tp) < 0) then
    289                
    290                 ! MG Mie approach - adjust density of sphere with D = D_characteristic to match particle density               
    291                 ! hp%rho_eff(tp,1:ns,iRe_type) = (6/pi)*hp%apm(tp)*(Di*1E-6)**(hp%bpm(tp)-3)   !MG Mie approach
    292                
    293                 ! as the particle size gets small it is possible that the mass to size relationship of
    294                 ! (given by power law in hclass.data) can produce impossible results
    295                 ! where the mass is larger than a solid sphere of ice. 
    296                 ! This loop ensures that no ice particle can have more mass/density larger than an ice sphere.
    297                 ! do i=1,ns
    298                 ! if(hp%rho_eff(tp,i,iRe_type) > 917 ) then
    299                 !       hp%rho_eff(tp,i,iRe_type) = 917
    300                 !endif
    301                 !enddo
    302 
    303                 ! alternative is to use equivalent volume spheres.
    304                 hp%rho_eff(tp,1:ns,iRe_type) = 917                              ! solid ice == equivalent volume approach
    305                 Deq = ( ( 6/pi*hp%apm(tp)/917 ) ** (1.0/3.0) ) * &
    306                            ( (Di*1E-6) ** (hp%bpm(tp)/3.0) )  * 1E6             ! Di now really Deq in microns.
    307                
     374      !     :: attenuation due to hydrometeors between radar and volume
     375      !
     376      ! NOTE old scheme integrates attenuation only for the layers ABOVE
     377      ! the current layer ... i.e. 1 to k-1 rather than 1 to k ...
     378      ! which may be a problem.   ROJ
     379      ! in the new scheme I assign half the attenuation to the current layer
     380      if(d_gate==1) then
     381        ! dheight calcuations assumes hgt_matrix points are the cell mid-points.
     382        if (k>2) then
     383          ! add to previous value to half of above layer + half of current layer
     384          a_to_vol(pr,k)=  a_to_vol(pr,k-1) + &
     385             (kr_vol(pr,k-1)+kr_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k))
     386        else
     387          a_to_vol(pr,k)=  kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1))
     388        endif
     389      else   ! d_gate==-1
     390        if(k<ngate) then
     391          ! add to previous value half of above layer + half of current layer
     392          a_to_vol(pr,k) = a_to_vol(pr,k+1) + &
     393              (kr_vol(pr,k+1)+kr_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k))
     394        else
     395          a_to_vol(pr,k)= kr_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1))
     396        endif
     397      endif
     398
     399      !     :: attenuation due to gaseous absorption between radar and volume
     400      if (g_to_vol_in_present) then
     401        g_to_vol(pr,k) = g_to_vol_in(pr,k)
     402      else
     403        if ( (hp%use_gas_abs == 1) .or. ((hp%use_gas_abs == 2) .and. (pr == 1)) ) then
     404          g_vol(pr,k) = gases(p_matrix(pr,k),t_kelvin,rh_matrix(pr,k),hp%freq)
     405          if (d_gate==1) then
     406            if (k>1) then
     407              ! add to previous value to half of above layer + half of current layer
     408              g_to_vol(pr,k) =  g_to_vol(pr,k-1) + &
     409                  0.5*(g_vol(pr,k-1)+g_vol(pr,k))*(hgt_matrix(pr,k-1)-hgt_matrix(pr,k))
    308410            else
    309 
    310                 ! hp%rho_eff(tp,1:ns,iRe_type) = hp%rho(tp)   !MG Mie approach
    311                
    312                 ! Equivalent volume sphere (solid ice rho_ice=917 kg/m^3).
    313                 hp%rho_eff(tp,1:ns,iRe_type) = 917
    314                 Deq=Di * ((hp%rho(tp)/917)**(1.0/3.0)) 
    315 
    316             endif
    317 
    318                 ! if using equivalent volume spheres
    319                 if (use_mie_table == 1) then
    320 
    321                         Ntemp=Ni
    322 
    323                         ! Find N(Di) from N(Deq) which we know
    324                         do i=1,ns
    325                                 j=infind(Deq,Di(i))
    326                                 Ni(i)=Ntemp(j)
    327                         enddo
    328                 else
    329                         ! just use Deq and D variable input to mie code
    330                         Di=Deq;
    331                 endif
    332 
    333           endif
    334           rhoi = hp%rho_eff(tp,1:ns,iRe_type)
    335          
    336 !         :: calculate effective reflectivity factor of volume
    337           if (use_mie_table == 1) then
    338          
    339             if ((hp%dtype(tp) == 4) .and. (hp%idd(tp) < 0)) then
    340               hp%idd(tp) = infind(mt%D,Di(1))
    341             endif
    342            
    343             if (phase == 0) then
    344            
    345               ! itt = infind(mt_ttl,t_matrix(pr,k))
    346               select case(hp%dtype(tp))
    347               case(4)
    348                 mt_qext(1) = mt%qext(hp%idd(tp),itt,1,iff)
    349                 mt_qbsca(1) = mt%qbsca(hp%idd(tp),itt,1,iff)
    350               case default
    351                 mt_qext = mt%qext(:,itt,1,iff)
    352                 mt_qbsca = mt%qbsca(:,itt,1,iff)
    353               end select
    354 
    355           call zeff(freq,Di,Ni,ns,k2,mt_ttl(itt),0,do_ray, &
    356                 ze,zr,kr,mt_qext,mt_qbsca,xx)
    357            
    358             else
    359 
    360               ! itt = infind(mt_tti,t_matrix(pr,k))
    361               select case(hp%dtype(tp))
    362               case(4)
    363                 if (hp%ifc(tp,1,iRe_type) < 0) then
    364                   hp%ifc(tp,1,iRe_type) = infind(mt%f,rhoi(1)/917.)
    365                 endif                 
    366                 mt_qext(1) = &
    367                   mt%qext(hp%idd(tp),itt+cnt_liq,hp%ifc(tp,1,iRe_type),iff)
    368                 mt_qbsca(1) = &
    369                   mt%qbsca(hp%idd(tp),itt+cnt_liq,hp%ifc(tp,1,iRe_type),iff)         
    370               case default
    371                 do i=1,ns
    372                   if (hp%ifc(tp,i,iRe_type) < 0) then
    373                     hp%ifc(tp,i,iRe_type) = infind(mt%f,rhoi(i)/917.)
    374                   endif       
    375                   mt_qext(i) = mt%qext(i,itt+cnt_liq,hp%ifc(tp,i,iRe_type),iff)
    376                   mt_qbsca(i) = mt%qbsca(i,itt+cnt_liq,hp%ifc(tp,i,iRe_type),iff)
    377                 enddo
    378               end select
    379 
    380                    call zeff(freq,Di,Ni,ns,k2,mt_tti(itt),1,do_ray, &
    381                 ze,zr,kr,mt_qext,mt_qbsca,xx)
    382 
    383             endif
    384 
    385           else
    386        
    387             xxa = -9.9
    388             call zeff(freq,Di,Ni,ns,k2,t_matrix(pr,k),phase,do_ray, &
    389               ze,zr,kr,xxa,xxa,rhoi)
    390 
    391              
    392           endif  ! end of use mie table
    393 
    394                 ! xxa = -9.9
    395                 !call zeff(freq,Di,Ni,ns,k2,t_matrix(pr,k),phase,do_ray, &
    396                 !       ze2,zr,kr2,xxa,xxa,rhoi)
    397 
    398                 ! if(abs(ze2-ze)/ze2 > 0.1) then
    399                 ! if(abs(kr2-kr)/kr2 > 0.1) then
    400                
    401                 ! write(*,*) pr,k,tp,ze2,ze2-ze,abs(ze2-ze)/ze2,itt+cnt_liq,iff
    402                 ! write(*,*) pr,k,tp,ze2,kr2,kr2-kr,abs(kr2-kr)/kr2
    403                 ! stop
    404 
    405                 !endif
    406 
    407           deallocate(Di,Ni,rhoi,xxa,Deq)
    408           if (use_mie_table == 1) deallocate(mt_qext,mt_qbsca,Ntemp)
    409 
    410           else ! can use z scaling
    411          
    412                 if( hp%dtype(tp)==2 .and. abs(hp%p2(tp)+1) < 1E-8 )  then
    413                  
    414                         ze = hp%Ze_scaled(tp,itt,iRe_type)
    415                         zr = hp%Zr_scaled(tp,itt,iRe_type)
    416                         kr = hp%kr_scaled(tp,itt,iRe_type)
    417 
    418                 else
    419                         scale_factor=rho_a*hm_matrix(tp,pr,k)
    420 
    421                         zr = hp%Zr_scaled(tp,itt,iRe_type) * scale_factor
    422                         ze = hp%Ze_scaled(tp,itt,iRe_type) * scale_factor
    423                         kr = hp%kr_scaled(tp,itt,iRe_type) * scale_factor       
    424                 endif
    425 
    426           endif  ! end z_scaling
    427  
    428           ! kr=0
    429 
    430           kr_vol(k) = kr_vol(k) + kr
    431           z_vol(k) = z_vol(k) + ze
    432           z_ray(k) = z_ray(k) + zr
    433        
    434           ! construct Ze_scaled, Zr_scaled, and kr_scaled ... if we can
    435           if( .not. hp%z_flag(tp,itt,iRe_type) .and. 1.eq.1 ) then
    436 
    437                 if( ( (hp%dtype(tp)==1 .or. hp%dtype(tp)==5 .or.  hp%dtype(tp)==2)  .and. abs(hp%p1(tp)+1) < 1E-8  ) .or. &
    438                     (  hp%dtype(tp)==3 .or. hp%dtype(tp)==4 )  &
    439                 ) then
    440 
    441                         scale_factor=rho_a*hm_matrix(tp,pr,k)
    442 
    443                         hp%Ze_scaled(tp,itt,iRe_type) = ze/ scale_factor
    444                         hp%Zr_scaled(tp,itt,iRe_type) = zr/ scale_factor
    445                         hp%kr_scaled(tp,itt,iRe_type) = kr/ scale_factor
    446 
    447                         hp%z_flag(tp,itt,iRe_type)=.True.
    448 
    449                 elseif( hp%dtype(tp)==2 .and. abs(hp%p2(tp)+1) < 1E-8 ) then
    450                  
    451                         hp%Ze_scaled(tp,itt,iRe_type) = ze
    452                         hp%Zr_scaled(tp,itt,iRe_type) = zr
    453                         hp%kr_scaled(tp,itt,iRe_type) = kr
    454 
    455                         hp%z_flag(tp,itt,iRe_type)=.True.
    456                 endif
    457 
    458           endif
    459 
    460         enddo   ! end loop of tp (hydrometeor type)
    461 
     411              g_to_vol(pr,k)=  0.5*g_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k+1))
     412            endif
     413          else   ! d_gate==-1
     414            if (k<ngate) then
     415              ! add to previous value to half of above layer + half of current layer
     416              g_to_vol(pr,k) = g_to_vol(pr,k+1) + &
     417                 0.5*(g_vol(pr,k+1)+g_vol(pr,k))*(hgt_matrix(pr,k+1)-hgt_matrix(pr,k))
     418            else
     419              g_to_vol(pr,k)= 0.5*g_vol(pr,k)*(hgt_matrix(pr,k)-hgt_matrix(pr,k-1))
     420            endif
     421          endif
     422        elseif(hp%use_gas_abs == 2) then
     423          ! using value calculated for the first column
     424          g_to_vol(pr,k) = g_to_vol(1,k)
     425        elseif (hp%use_gas_abs == 0) then
     426          g_to_vol(pr,k) = 0
     427        endif
     428      endif
     429
     430      ! Compute Rayleigh reflectivity, and full, attenuated reflectivity
     431      if ((hp%do_ray == 1) .and. (z_ray(pr,k) > 0)) then
     432        Ze_ray(pr,k) = 10*log10(z_ray(pr,k))
    462433      else
    463 !     :: volume is hydrometeor-free
    464        
    465         kr_vol(k) = 0
    466         z_vol(k) = -999
    467         z_ray(k) = -999
    468        
     434        Ze_ray(pr,k) = undef
    469435      endif
    470 
    471 !     :: attenuation due to hydrometeors between radar and volume
    472       a_to_vol(k) = 2*path_integral(kr_vol,hgt_matrix(pr,:),1,k-1)
    473      
    474 !     :: attenuation due to gaseous absorption between radar and volume
    475       if (g_to_vol_in_present) then
    476         g_to_vol(k) = g_to_vol_in(k,pr)
     436      if (z_vol(pr,k) > 0) then
     437        Ze_non(pr,k) = 10*log10(z_vol(pr,k))
     438        dBZe(pr,k) = Ze_non(pr,k)-a_to_vol(pr,k)-g_to_vol(pr,k)
    477439      else
    478         if ( (use_gas_abs == 1) .or. ((use_gas_abs == 2) .and. (pr == 1)) )  then
    479             g_vol(k) = gases(p_matrix(pr,k),t_matrix(pr,k)+273.15, &
    480             rh_matrix(pr,k),freq)
    481             g_to_vol(k) = path_integral(g_vol,hgt_matrix(pr,:),1,k-1)
    482         elseif (use_gas_abs == 0) then
    483             g_to_vol(k) = 0
    484         endif 
     440        dBZe(pr,k) = undef
     441        Ze_non(pr,k) = undef
    485442      endif
    486    
    487 !      kr_matrix(pr,:)=kr_vol
    488 
    489 !     :: store results in matrix for return to calling program
    490       h_atten_to_vol(pr,k)=a_to_vol(k)
    491       g_atten_to_vol(pr,k)=g_to_vol(k)
    492       if ((do_ray == 1) .and. (z_ray(k) > 0)) then
    493         Ze_ray(pr,k) = 10*log10(z_ray(k))
    494       else
    495         Ze_ray(pr,k) = -999
    496       endif
    497       if (z_vol(k) > 0) then
    498         dBZe(pr,k) = 10*log10(z_vol(k))-a_to_vol(k)-g_to_vol(k)
    499         Ze_non(pr,k) = 10*log10(z_vol(k))
    500       else
    501         dBZe(pr,k) = -999
    502         Ze_non(pr,k) = -999
    503       endif
    504      
    505     enddo       ! end loop of k (range gate)
    506     ! Output array with gaseous absorption
    507     if (g_to_vol_out_present) g_to_vol_out(:,pr) = g_to_vol
    508   enddo         ! end loop over pr (profile) 
     443
     444    enddo   ! end loop over pr (profile)
     445
     446  enddo ! end loop of k (range gate)
     447
     448  ! Output array with gaseous absorption
     449  if (g_to_vol_out_present) g_to_vol_out = g_to_vol
     450
     451  ! save any updates made
     452  if (hp%update_scale_LUTs) call save_scale_LUTs(hp)
    509453
    510454  end subroutine radar_simulator
    511  
Note: See TracChangeset for help on using the changeset viewer.