phys/module_mp_SBM_polar_radar.F

   1 #if( BUILD_SBM_FAST != 1)
   2       MODULE module_mp_SBM_polar_radar
   3       CONTAINS
   4       SUBROUTINE SBM_polar_radar
   5          REAL :: dummy
   6          dummy = 1
   7       END SUBROUTINE SBM_polar_radar
   8       END MODULE module_mp_SBM_polar_radar
   9 #else
  10 !******************
  11 module scatt_tables
  12 ! JCS - This module pertains to the reading of scattering amplitude files
  13
  14 !use microprm
  15
  16 implicit none
  17
  18 private
  19 public :: faf1,fbf1,fab1,fbb1,         &
  20         ! faf1fd,fbf1fd,fab1fd,fbb1fd, &
  21         ! faf2d,fbf2d,fab2d,fbb2d,     &
  22         ! faf2p,fbf2p,fab2p,fbb2p,     &
  23         ! faf2c,fbf2c,fab2c,fbb2c,     &
  24           faf3,fbf3,fab3,fbb3,         &
  25           faf4,fbf4,fab4,fbb4,         &
  26           faf5,fbf5,fab5,fbb5,         &
  27           LOAD_TABLES,                 &
  28           temps_water,temps_fd,temps_crystals,      &
  29           temps_snow,temps_graupel,temps_hail,      &
  30           fws_fd,fws_crystals,fws_snow,             &
  31           fws_graupel,fws_hail,                     &
  32           usetables,                                &
  33           twolayer_hail,twolayer_graupel,twolayer_fd,twolayer_snow,rpquada,usequad
  34
  35 SAVE ! [KS >> This "SAVE" possibly interfering]
  36
  37 ! JCS -- if usetables is TRUE, this module will read the precomputed scattering amplitudes.
  38 ! If usetables is 0, this module will do nothing, and the program will
  39 ! calculate the scattering amplitudes as necessary within the program (that is,
  40 ! no lookup tables will be used). If usetables is 1, we'll use precomputed
  41 ! scattering amplitudes. usetables(water,fd,crystals,snow,graupel,hail)
  42 integer, dimension(6) :: usetables = (/1,0,0,1,1,1/)
  43 ! JCS -- If set to 1, the two-layer T-matrix scattering code will be used where
  44 ! necessary. If 0, then we'll only use the homogeneous-mixture T-matrix code.
  45 integer :: twolayer_hail = 1
  46 integer :: twolayer_graupel = 1
  47 integer :: twolayer_fd = 1
  48 integer :: twolayer_snow = 1
  49
  50 ! JCS - Use quad precision for two-layer calculations for large sizes?
  51 logical,parameter :: usequad = .true.
  52 ! JCS - If usequad is true, then rpquada will be used to define the rp at which
  53 ! quad-precision 2-layer t-matrix will be called (that is, quad-precision is
  54 ! used for rp >= rpquada). rpquada(water, fd, snow/crystals, gruapel, hail)
  55 double precision, dimension(5) :: rpquada = (/10.0d1,2.0d0,1.5d0,2.0d0,2.0d0/)
  56
  57 ! JCS -- in the current version of HUCM, each species has the same number of
  58 ! bins (i.e., NKR, which is set in the microprm.F90 file/module).
  59 ! >> KS (comment-out) : integer, parameter :: nbins=NKR
  60
  61 ! JCS -- each hydrometeor species will have a 3-dimensional table sized
  62 ! NKR x ntemps x nfws (that is, number of bins by number of temperatures by
  63 ! number of water fractions).
  64 ! JCS -- arrays are ordered as (water,fd,snow/crystals,graupel,hail)
  65 integer, dimension(5),parameter :: tstart = (/-20,-20,-20,-20,-20/), ntemps = (/61,31,31,61,61/),  &
  66                                    dtemp = (/1,1,1,1,1/), nfws = (/1,101,101,101,101/)
  67
  68 ! >> [KS] integer, allocatable :: temps_water(:), temps_fd(:), temps_crystals(:), temps_snow(:), temps_graupel(:), temps_hail(:)
  69 integer :: i,ios,iiwl,ispecies
  70 integer, parameter, dimension(ntemps(1)) :: temps_water=(/ (dtemp(1)*(i-1)+tstart(1),i=1,ntemps(1) )/)
  71 integer, parameter, dimension(ntemps(2)) :: temps_fd=(/(dtemp(2)*(i-1)+tstart(2),i=1,ntemps(2))/)
  72 integer, parameter, dimension(ntemps(3)) :: temps_crystals=(/(dtemp(3)*(i-1)+tstart(3),i=1,ntemps(3))/)
  73 integer, parameter, dimension(ntemps(3)) :: temps_snow=(/(dtemp(3)*(i-1)+tstart(3),i=1,ntemps(3))/)
  74 integer, parameter, dimension(ntemps(4)) :: temps_graupel=(/(dtemp(4)*(i-1)+tstart(4),i=1,ntemps(4))/)
  75 integer, parameter, dimension(ntemps(5)) :: temps_hail=(/(dtemp(5)*(i-1)+tstart(5),i=1,ntemps(5))/)
  76
  77 ! JCS - If fvw=1.0 (i.e., 100%), then it should be considered rain.
  78 ! Units are decimal fraction from 0.0 to 1.0
  79 real :: fws_water=1.0
  80 ! >> [KS] real, allocatable :: fws_fd(:), fws_crystals(:), fws_snow(:), fws_graupel(:), fws_hail(:)
  81 real, parameter, dimension(nfws(2)) :: fws_fd=(/(1.0/(nfws(2)-1)*(i-1),i=1,nfws(2))/)
  82 real, parameter, dimension(nfws(3)) :: fws_crystals=(/(1.0/(nfws(3)-1)*(i-1),i=1,nfws(3))/)
  83 real, parameter, dimension(nfws(3)) :: fws_snow=(/(1.0/(nfws(3)-1)*(i-1),i=1,nfws(3))/)
  84 real, parameter, dimension(nfws(4)) :: fws_graupel=(/(1.0/(nfws(4)-1)*(i-1),i=1,nfws(4))/)
  85 real, parameter, dimension(nfws(5)) :: fws_hail=(/(1.0/(nfws(5)-1)*(i-1),i=1,nfws(5))/)
  86
  87 ! JCS - Array of wavelengths used in the polarimetric emulator/polar_hucm.F90
  88 ! Wavelengths should be in units of cm. The number of wavelengths must match the
  89 ! number of filesnames and the size of the FILENAMES array
  90 INTEGER,parameter :: nwavelengths = 1
  91 DOUBLE PRECISION :: WAVELENGTH1, WAVELENGTH2, WAVELENGTH3
  92 DOUBLE PRECISION, DIMENSION(3),parameter :: WAVELENGTHS = (/11.0D0, 5.5D0, 3.2D0/)
  93 CHARACTER(LEN=20),parameter :: OUTFILENAME1='GRADS_MOV_SBAND', OUTFILENAME2='GRADS_MOV_CBAND', &
  94                                 OUTFILENAME3='GRADS_MOV_XBAND'
  95
  96 ! >> [KS] CHARACTER(LEN=20), ALLOCATABLE :: FILENAMES(:)
  97 CHARACTER(LEN=20),parameter,dimension(nwavelengths) :: FILENAMES=(/OUTFILENAME1/)
  98
  99 CHARACTER(LEN=256),parameter :: scattering_dir_prefix = 'scattering_tables_2layer_high_quad_1dT_1%fw'
 100 ! >> [KS] CHARACTER(LEN=256), ALLOCATABLE :: scattering_dir(:)
 101 CHARACTER(LEN=256),dimension(nwavelengths) :: scattering_dir
 102
 103 CHARACTER(Len=3) :: wlstr
 104
 105 ! JCS - below are the tables that will hold the scattering amplitudes
 106 ! They are named as f(a for horizontal, b for vertical)(f for forward, b for
 107 ! backward),(1 for water, 1fd for freezing drops, 2d for dendrites, 2p for plates,
 108 ! 2c for columns, 3 for snow aggregates, 4 for graupel, and 5 for hail)
 109 double complex, allocatable :: faf1(:,:,:,:),fbf1(:,:,:,:),fab1(:,:,:,:),fbb1(:,:,:,:)
 110 double complex, allocatable :: faf1fd(:,:,:,:),fbf1fd(:,:,:,:),fab1fd(:,:,:,:),fbb1fd(:,:,:,:)
 111 double complex, allocatable :: faf2d(:,:,:,:),fbf2d(:,:,:,:),fab2d(:,:,:,:),fbb2d(:,:,:,:)
 112 double complex, allocatable :: faf2p(:,:,:,:),fbf2p(:,:,:,:),fab2p(:,:,:,:),fbb2p(:,:,:,:)
 113 double complex, allocatable :: faf2c(:,:,:,:),fbf2c(:,:,:,:),fab2c(:,:,:,:),fbb2c(:,:,:,:)
 114 double complex, allocatable :: faf3(:,:,:,:),fbf3(:,:,:,:),fab3(:,:,:,:),fbb3(:,:,:,:)
 115 double complex, allocatable :: faf4(:,:,:,:),fbf4(:,:,:,:),fab4(:,:,:,:),fbb4(:,:,:,:)
 116 double complex, allocatable :: faf5(:,:,:,:),fbf5(:,:,:,:),fab5(:,:,:,:),fbb5(:,:,:,:)
 117
 118 integer, dimension(1) :: itemp, infw
 119 ! >> [KS] integer :: ispecies, i, ios, iiwl
 120
 121 !NAMELIST /scatttables/ usetables,twolayer_hail,twolayer_graupel,twolayer_fd,twolayer_snow, &
 122 !                       usequad,rpquada,tstart,ntemps,dtemp,nfws, &
 123 !                       nwavelengths, wavelengths, outfilename1, outfilename2, outfilename3, &
 124 !                       scattering_dir_prefix
 125
 126 CONTAINS
 127
 128 SUBROUTINE LOAD_TABLES(nbins)
 129
 130 implicit none
 131
 132 integer istatus
 133 character*256 :: header,header2
 134 character*256 :: fname
 135 integer :: i,j,k
 136 character*3 :: temp, fw
 137 integer,intent(in) :: nbins   ! >> (KS)
 138 real,dimension(nbins) :: m    ! >> (KS)
 139
 140
 141 !OPEN(101,FILE='scatt_tables.input',STATUS="old")
 142 !    read(101,scatttables)
 143 !CLOSE(101)
 144 !print *,wavelengths
 145
 146 !>>ALLOCATE(temps_water(ntemps(1)),stat=istatus)
 147 !>>ALLOCATE(temps_fd(ntemps(2)),stat=istatus)
 148 !>>ALLOCATE(temps_crystals(ntemps(3)),stat=istatus)
 149 !>>ALLOCATE(temps_snow(ntemps(3)),stat=istatus)
 150 !>>ALLOCATE(temps_graupel(ntemps(4)),stat=istatus)
 151 !>>ALLOCATE(temps_hail(ntemps(5)),stat=istatus)
 152
 153 !>>ALLOCATE(fws_fd(nfws(2)),stat=istatus)
 154 !>>ALLOCATE(fws_crystals(nfws(3)),stat=istatus)
 155 !>>ALLOCATE(fws_snow(nfws(3)),stat=istatus)
 156 !>>ALLOCATE(fws_graupel(nfws(4)),stat=istatus)
 157 !>>ALLOCATE(fws_hail(nfws(5)),stat=istatus)
 158
 159 !>>temps_water=(/ (dtemp(1)*(i-1)+tstart(1),i=1,ntemps(1) )/)
 160 !>>temps_fd=(/(dtemp(2)*(i-1)+tstart(2),i=1,ntemps(2))/)
 161 !>>temps_crystals=(/(dtemp(3)*(i-1)+tstart(3),i=1,ntemps(3))/)
 162 !>>temps_snow=(/(dtemp(3)*(i-1)+tstart(3),i=1,ntemps(3))/)
 163 !>>temps_graupel=(/(dtemp(4)*(i-1)+tstart(4),i=1,ntemps(4))/)
 164 !>>temps_hail=(/(dtemp(5)*(i-1)+tstart(5),i=1,ntemps(5))/)
 165
 166 !>>fws_fd=(/(1.0/(nfws(2)-1)*(i-1),i=1,nfws(2))/)
 167 !>>fws_crystals=(/(1.0/(nfws(3)-1)*(i-1),i=1,nfws(3))/)
 168 !>>fws_snow=(/(1.0/(nfws(3)-1)*(i-1),i=1,nfws(3))/)
 169 !>>fws_graupel=(/(1.0/(nfws(4)-1)*(i-1),i=1,nfws(4))/)
 170 !>>fws_hail=(/(1.0/(nfws(5)-1)*(i-1),i=1,nfws(5))/)
 171
 172 !>>ALLOCATE(FILENAMES(nwavelengths),stat=istatus)
 173 !ALLOCATE(WAVELENGTHS(nwavelengths),stat=istatus)
 174 !>>ALLOCATE(scattering_dir(nwavelengths),stat=istatus)
 175
 176 !>>FILENAMES=(/OUTFILENAME1,OUTFILENAME2,OUTFILENAME3/)
 177 !WAVELENGTHS=(/WAVELENGTH1,WAVELENGTH2,WAVELENGTH3/)
 178
 179 do iiwl=1,nwavelengths
 180   write(wlstr,'(I3.3)') int(WAVELENGTHS(iiwl)*10.0d0)
 181   scattering_dir(iiwl)=TRIM(scattering_dir_prefix)//'_'//wlstr//'/'
 182   WRITE(*,*) 'scattering input directory is source/',TRIM(scattering_dir(iiwl))
 183 enddo
 184
 185 DO ispecies=1,size(usetables)
 186   if((ispecies==1) .AND. (usetables(ispecies)==1) ) then ! rain
 187       WRITE(*,*) 'READING SCATTERING TABLES: RAIN'
 188       ALLOCATE(faf1(nbins,nfws(1),ntemps(1),nwavelengths),stat=istatus)
 189       ALLOCATE(fbf1(nbins,nfws(1),ntemps(1),nwavelengths),stat=istatus)
 190       ALLOCATE(fab1(nbins,nfws(1),ntemps(1),nwavelengths),stat=istatus)
 191       ALLOCATE(fbb1(nbins,nfws(1),ntemps(1),nwavelengths),stat=istatus)
 192       do iiwl=1,nwavelengths
 193         do k=1,ntemps(1)
 194           write(temp,"(SP,I3.2)") temps_water(k)
 195           !>>fname='source/'//TRIM(scattering_dir(iiwl))//'/RAIN_'//temp//'C_100fvw.sct'
 196           fname=TRIM(scattering_dir(iiwl))//'/RAIN_'//temp//'C_100fvw.sct'
 197 !!          WRITE(*,*) TRIM(fname)
 198           open(unit=1,file=fname,status="old",form="formatted",iostat=ios)
 199           if(ios.eq.0) then
 200               read(1,*) header
 201               read(1,*) header2
 202               do i=1,nbins
 203                 read(1,'(9e13.5)') m(i), fab1(i,1,k,iiwl), fbb1(i,1,k,iiwl), &
 204                                         faf1(i,1,k,iiwl), fbf1(i,1,k,iiwl)
 205               enddo
 206               close(1)
 207           else
 208               WRITE(*,*) '*****PROBLEM READING RAIN SCATTERING AMPLITUDE FILE*****'
 209               WRITE(*,*) TRIM(fname)
 210 !!              WRITE(*,*) 'Temp=',TRIM(temp),' C, fvw=100'
 211               WRITE(*,*) '*****NO LOOKUP TABLES FOR RAIN WILL BE USED*****'
 212               usetables(1)=0
 213           endif
 214         enddo
 215       enddo
 216   elseif(ispecies==2 .AND. usetables(ispecies)==1) then ! fd
 217       WRITE(*,*) 'READING SCATTERING TABLES: FD'
 218       ALLOCATE(faf1fd(nbins,nfws(2),ntemps(2),nwavelengths),stat=istatus)
 219       ALLOCATE(fbf1fd(nbins,nfws(2),ntemps(2),nwavelengths),stat=istatus)
 220       ALLOCATE(fab1fd(nbins,nfws(2),ntemps(2),nwavelengths),stat=istatus)
 221       ALLOCATE(fbb1fd(nbins,nfws(2),ntemps(2),nwavelengths),stat=istatus)
 222       do iiwl=1,nwavelengths
 223        do k=1,ntemps(2)
 224         write(temp,"(SP,I3.2)") temps_fd(k)
 225         do j=1,nfws(2)
 226           write(fw,"(I3.3)") NINT(fws_fd(j)*100)
 227           !>>fname='source/'//TRIM(scattering_dir(iiwl))//'/FD_'//temp//'C_'//fw//'fvw.sct'
 228           fname=TRIM(scattering_dir(iiwl))//'/FD_'//temp//'C_'//fw//'fvw.sct'
 229           open(unit=1,file=fname,status="old",form="formatted",iostat=ios)
 230           if(ios.eq.0) then
 231               read(1,*) header
 232               read(1,*) header2
 233               do i=1,nbins
 234                 read(1,'(9e13.5)') m(i), fab1fd(i,j,k,iiwl), fbb1fd(i,j,k,iiwl), &
 235                     faf1fd(i,j,k,iiwl), fbf1fd(i,j,k,iiwl)
 236               enddo
 237               close(1)
 238           else
 239               WRITE(*,*) '*****PROBLEM READING FD SCATTERING AMPLITUDE FILE*****'
 240               WRITE(*,*) TRIM(fname)
 241 !!              WRITE(*,*) 'Temp=',TRIM(temp),' C, fvw=',TRIM(fw)
 242               WRITE(*,*) '*****NO LOOKUP TABLES FOR FD WILL BE USED*****'
 243               usetables(2)=0
 244           endif
 245         enddo
 246        enddo
 247       enddo
 248   elseif(ispecies==3 .AND. usetables(ispecies)==1) then ! ice crystals (plates, dendrites, columns)
 249       WRITE(*,*) 'READING SCATTERING TABLES: ICE CRYSTALS'
 250       ALLOCATE(faf2d(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 251       ALLOCATE(fbf2d(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 252       ALLOCATE(fab2d(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 253       ALLOCATE(fbb2d(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 254       ALLOCATE(faf2p(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 255       ALLOCATE(fbf2p(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 256       ALLOCATE(fab2p(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 257       ALLOCATE(fbb2p(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 258       ALLOCATE(faf2c(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 259       ALLOCATE(fbf2c(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 260       ALLOCATE(fab2c(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 261       ALLOCATE(fbb2c(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 262       do iiwl=1,nwavelengths
 263        do k=1,ntemps(3)
 264         write(temp,"(SP,I3.2)") temps_crystals(k)
 265         do j=1,nfws(3)
 266           write(fw,"(I3.3)") NINT(fws_crystals(j)*100)
 267           !>>fname='source/'//TRIM(scattering_dir(iiwl))//'/DENDRITES_'//temp//'C_'//fw//'fvw.sct'
 268           fname=TRIM(scattering_dir(iiwl))//'/DENDRITES_'//temp//'C_'//fw//'fvw.sct'
 269           open(unit=1,file=fname,status="old",form="formatted",iostat=ios)
 270           if(ios.eq.0) then
 271               read(1,*) header
 272               read(1,*) header2
 273               do i=1,nbins
 274                 read(1,'(9e13.5)') m(i), fab2d(i,j,k,iiwl), fbb2d(i,j,k,iiwl), &
 275                     faf2d(i,j,k,iiwl), fbf2d(i,j,k,iiwl)
 276               enddo
 277               close(1)
 278           else
 279               WRITE(*,*) '*****PROBLEM READING DENDRITES SCATTERING AMPLITUDE FILE*****'
 280               WRITE(*,*) 'Temp=',TRIM(temp),' C, fvw=',TRIM(fw)
 281               WRITE(*,*) '*****NO LOOKUP TABLES FOR DENDRITES WILL BE USED*****'
 282               usetables(3)=0
 283           endif
 284         enddo
 285        enddo
 286        do k=1,ntemps(3)
 287         write(temp,"(SP,I3.2)") temps_crystals(k)
 288         do j=1,nfws(3)
 289           write(fw,"(I3.3)") NINT(fws_crystals(j)*100)
 290           !>>fname='source/'//TRIM(scattering_dir(iiwl))//'/PLATES_'//temp//'C_'//fw//'fvw.sct'
 291           fname=TRIM(scattering_dir(iiwl))//'/PLATES_'//temp//'C_'//fw//'fvw.sct'
 292           open(unit=1,file=fname,status="old",form="formatted",iostat=ios)
 293           if(ios.eq.0) then
 294               read(1,*) header
 295               read(1,*) header2
 296               do i=1,nbins
 297                 read(1,'(f5.2,8e13.5)') m(i), fab2p(i,j,k,iiwl), fbb2p(i,j,k,iiwl), &
 298                     faf2p(i,j,k,iiwl), fbf2p(i,j,k,iiwl)
 299               enddo
 300               close(1)
 301           else
 302               WRITE(*,*) '*****PROBLEM READING PLATES SCATTERING AMPLITUDE FILE*****'
 303               WRITE(*,*) 'Temp=',TRIM(temp),' C, fvw=',TRIM(fw)
 304               WRITE(*,*) '*****NO LOOKUP TABLES FOR PLATES WILL BE USED*****'
 305               usetables(3)=0
 306           endif
 307         enddo
 308        enddo
 309        do k=1,ntemps(3)
 310         write(temp,"(SP,I3.2)") temps_crystals(k)
 311         do j=1,nfws(3)
 312           write(fw,"(I3.3)") NINT(fws_crystals(j)*100)
 313           !fname='source/'//TRIM(scattering_dir(iiwl))//'/COLUMNS_'//temp//'C_'//fw//'fvw.sct'
 314           fname=TRIM(scattering_dir(iiwl))//'/COLUMNS_'//temp//'C_'//fw//'fvw.sct'
 315           open(unit=1,file=fname,status="old",form="formatted",iostat=ios)
 316           if(ios.eq.0) then
 317               read(1,*) header
 318               read(1,*) header2
 319               do i=1,nbins
 320                 read(1,'(9e13.5)') m(i), fab2c(i,j,k,iiwl), fbb2c(i,j,k,iiwl), &
 321                     faf2c(i,j,k,iiwl), fbf2c(i,j,k,iiwl)
 322               enddo
 323               close(1)
 324           else
 325               WRITE(*,*) '*****PROBLEM READING COLUMNS SCATTERING AMPLITUDE FILE*****'
 326               WRITE(*,*) 'Temp=',TRIM(temp),' C, fvw=',TRIM(fw)
 327               WRITE(*,*) '*****NO LOOKUP TABLES FOR COLUMNS WILL BE USED*****'
 328               usetables(3)=0
 329           endif
 330         enddo
 331        enddo
 332       enddo
 333   elseif(ispecies==4 .AND. usetables(ispecies)==1) then ! snow (aggregates)
 334       WRITE(*,*) 'READING SCATTERING TABLES: SNOW'
 335       ALLOCATE(faf3(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 336       ALLOCATE(fbf3(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 337       ALLOCATE(fab3(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 338       ALLOCATE(fbb3(nbins,nfws(3),ntemps(3),nwavelengths),stat=istatus)
 339       do iiwl=1,nwavelengths
 340        do k=1,ntemps(3)
 341         write(temp,"(SP,I3.2)") temps_snow(k)
 342         do j=1,nfws(3)
 343           write(fw,"(I3.3)") NINT(fws_snow(j)*100)
 344           !>>fname='source/'//TRIM(scattering_dir(iiwl))//'/SNOW_'//temp//'C_'//fw//'fvw.sct'
 345           fname=TRIM(scattering_dir(iiwl))//'/SNOW_'//temp//'C_'//fw//'fvw.sct'
 346           open(unit=1,file=fname,status="old",form="formatted",iostat=ios)
 347           if(ios.eq.0) then
 348               read(1,*) header
 349               read(1,*) header2
 350               do i=1,nbins
 351                 read(1,'(9e13.5)') m(i), fab3(i,j,k,iiwl), fbb3(i,j,k,iiwl), &
 352                     faf3(i,j,k,iiwl), fbf3(i,j,k,iiwl)
 353               enddo
 354               close(1)
 355           else
 356               WRITE(*,*) '*****PROBLEM READING SNOW SCATTERING AMPLITUDE FILE*****'
 357               WRITE(*,*) TRIM(fname)
 358 !!              WRITE(*,*) 'Temp=',TRIM(temp),' C, fvw=',TRIM(fw)
 359               WRITE(*,*) '*****NO LOOKUP TABLES FOR SNOW WILL BE USED*****'
 360               usetables(4)=0
 361           endif
 362         enddo
 363        enddo
 364       enddo
 365   elseif(ispecies==5 .AND. usetables(ispecies)==1) then ! graupel
 366       WRITE(*,*) 'READING SCATTERING TABLES: GRAUPEL'
 367       ALLOCATE(faf4(nbins,nfws(4),ntemps(4),nwavelengths),stat=istatus)
 368       ALLOCATE(fbf4(nbins,nfws(4),ntemps(4),nwavelengths),stat=istatus)
 369       ALLOCATE(fab4(nbins,nfws(4),ntemps(4),nwavelengths),stat=istatus)
 370       ALLOCATE(fbb4(nbins,nfws(4),ntemps(4),nwavelengths),stat=istatus)
 371       do iiwl=1,nwavelengths
 372        do k=1,ntemps(4)
 373         write(temp,"(SP,I3.2)") temps_graupel(k)
 374         do j=1,nfws(4)
 375           write(fw,"(I3.3)") NINT(fws_graupel(j)*100)
 376           !>>fname='source/'//TRIM(scattering_dir(iiwl))//'/GRAUPEL_'//temp//'C_'//fw//'fvw.sct'
 377           fname=TRIM(scattering_dir(iiwl))//'/GRAUPEL_'//temp//'C_'//fw//'fvw.sct'
 378           open(unit=1,file=fname,status="old",form="formatted",iostat=ios)
 379           if(ios.eq.0) then
 380               read(1,*) header
 381               read(1,*) header2
 382               do i=1,nbins
 383                 read(1,'(9e13.5)') m(i), fab4(i,j,k,iiwl), fbb4(i,j,k,iiwl), &
 384                     faf4(i,j,k,iiwl), fbf4(i,j,k,iiwl)
 385               enddo
 386               close(1)
 387           else
 388               WRITE(*,*) '*****PROBLEM READING GRAUPEL SCATTERING AMPLITUDE FILE*****'
 389               WRITE(*,*) TRIM(fname)
 390 !!              WRITE(*,*) 'Temp=',TRIM(temp),' C, fvw=',TRIM(fw)
 391               WRITE(*,*) '*****NO LOOKUP TABLES FOR GRAUPEL WILL BE USED*****'
 392               usetables(5)=0
 393           endif
 394         enddo
 395        enddo
 396       enddo
 397   elseif(ispecies==6 .AND. usetables(ispecies)==1) then ! hail
 398       WRITE(*,*) 'READING SCATTERING TABLES: HAIL'
 399       ALLOCATE(faf5(nbins,nfws(5),ntemps(5),nwavelengths),stat=istatus)
 400       ALLOCATE(fbf5(nbins,nfws(5),ntemps(5),nwavelengths),stat=istatus)
 401       ALLOCATE(fab5(nbins,nfws(5),ntemps(5),nwavelengths),stat=istatus)
 402       ALLOCATE(fbb5(nbins,nfws(5),ntemps(5),nwavelengths),stat=istatus)
 403       do iiwl=1,nwavelengths
 404        do k=1,ntemps(5)
 405         write(temp,"(SP,I3.2)") temps_hail(k)
 406         do j=1,nfws(5)
 407           write(fw,"(I3.3)") NINT(fws_hail(j)*100)
 408           !>>fname='source/'//TRIM(scattering_dir(iiwl))//'/HAIL_'//temp//'C_'//fw//'fvw.sct'
 409           fname=TRIM(scattering_dir(iiwl))//'/HAIL_'//temp//'C_'//fw//'fvw.sct'
 410           open(unit=1,file=fname,status="old",form="formatted",iostat=ios)
 411           if(ios.eq.0) then
 412               read(1,*) header
 413               read(1,*) header2
 414               do i=1,nbins
 415                 read(1,'(9e13.5)') m(i), fab5(i,j,k,iiwl), fbb5(i,j,k,iiwl), &
 416                     faf5(i,j,k,iiwl), fbf5(i,j,k,iiwl)
 417               enddo
 418               close(1)
 419           else
 420               WRITE(*,*) '*****PROBLEM READING HAIL SCATTERING AMPLITUDE FILE*****'
 421               WRITE(*,*) TRIM(fname)
 422 !!              WRITE(*,*) 'Temp=',TRIM(temp),' C, fvw=',TRIM(fw)
 423               WRITE(*,*) '*****NO LOOKUP TABLES FOR HAIL WILL BE USED*****'
 424               usetables(6)=0
 425           endif
 426         enddo
 427        enddo
 428       enddo
 429   endif
 430 enddo
 431
 432 END SUBROUTINE LOAD_TABLES
 433
 434 SUBROUTINE CHECK_ALLOCATION_STATUS(istatus)
 435 implicit none
 436 integer :: istatus
 437
 438 END SUBROUTINE CHECK_ALLOCATION_STATUS
 439
 440
 441 END MODULE scatt_tables
 442 ! +----------------------------------------------------------------------------+
 443 ! +----------------------------------------------------------------------------+
 444 MODULE module_mp_SBM_polar_radar
 445
 446 ! Kind paramater
 447 INTEGER, PARAMETER, PRIVATE:: R8SIZE = 8
 448 INTEGER, PARAMETER, PRIVATE:: R4SIZE = 4
 449
 450 private
 451 public polar_hucm
 452
 453      !  Parameter (NPN1=100, NPNG1=500, NPNG2=2*NPNG1, NPN2=2*NPN1,NPL=NPN2+1, NPN3=NPN1+1,  NPN4=NPN1, NPN5=2*NPN4, NPN6=NPN4+1)
 454       LOGICAL, PRIVATE,PARAMETER :: TRANSMISSION=.FALSE.
 455      ! INTEGER, PRIVATE,PARAMETER :: ICEMAX = 3, NKR_43Bins = 43, NKR_33Bins = 33
 456      ! INTEGER, PRIVATE,PARAMETER :: ICEMAX=3,NKR=33
 457      ! INTEGER :: NKR
 458        CONTAINS
 459
 460
 461        subroutine polar_hucm &
 462                                   (FF1,FF2,FF3,FF4,FF5,FF1_FD,                                                          &
 463                                    FL3,FL4,FL5,FL1_FD,                                                                          &
 464                                    bulk,temp,RORD,wavelength,iwl,distance,                                      &
 465                                    dx,dy,zmks_1d,                                                                                       &
 466                                    out1,out2,out3,out4,out5,out6,out7,out8,out9,                        &
 467                                    bin_mass,tab_colum,tab_dendr, tab_snow, bin_log,             &
 468                                    ijk,kx,ky,kz,kts,kte,number_bin,ICEMAX,icloud,itimestep, &
 469                                    faf1,fbf1,fab1,fbb1,                                                                         &
 470                                    !faf1fd,fbf1fd,fab1fd,fbb1fd,                                                        &
 471                                    !faf2d,fbf2d,fab2d,fbb2d,                                                            &
 472                                    !faf2p,fbf2p,fab2p,fbb2p,                                                            &
 473                                    !faf2c,fbf2c,fab2c,fbb2c,                                                            &
 474                                         faf3,fbf3,fab3,fbb3,                                                            &
 475                                         faf4,fbf4,fab4,fbb4,                                                            &
 476                                         faf5,fbf5,fab5,fbb5,                                                            &
 477                                         temps_water,temps_fd,temps_crystals,                                    &
 478                                         temps_snow,temps_graupel,temps_hail,                                    &
 479                                         fws_fd,fws_crystals,fws_snow,                                                   &
 480                                         fws_graupel,fws_hail,usetables)
 481
 482 !**** *****************************************
 483 !     temperature          Celsius degree
 484 !     wavelength           cm
 485 !     density              g/cm^3
 486 !     equivolume diameter  mm
 487 !     amplitudes           mm
 488 !     mass                 g
 489 !     dx                   m
 490 !     dz                   m
 491 !     distance             m
 492 !     elevation            degree
 493 ! &
 494 !**** ******************************************
 495
 496           implicit none
 497
 498 ! ### (KS) : Interface Vars.
 499
 500           integer,intent(in) :: number_bin, icemax, kte, kts, kz, ky, kx, ijk, icloud, itimestep, iwl
 501           real(kind=r8size),intent(in) :: zmks_1d(KTS:KTE), bin_mass(number_bin), tab_colum(number_bin), tab_dendr(number_bin), &
 502                                                      tab_snow(number_bin),bin_log(number_bin), bulk(number_bin), temp, RORD, wavelength,   &
 503                                                      distance, dx,dy
 504           real(kind=r8size),intent(out) ::  out1(10), out2(10), out3(10), out4(10), &
 505                                                                                            out5(10),out6(10), out7(10),out8(10), out9(10)
 506           real(kind=r8size),intent(inout) :: FF1(number_bin),FF2(number_bin,ICEMAX),FF3(number_bin),FF4(number_bin),FF5(number_bin), &
 507                                                                           FF1_FD(number_bin), FL3(number_bin),FL4(number_bin),FL5(number_bin), &
 508                                                   FL1_FD(number_bin)
 509           double complex,intent(in), dimension(:,:,:,:) :: faf1,fbf1,fab1,fbb1,                                                                         &
 510                                                                                                           !faf1fd,fbf1fd,fab1fd,fbb1fd,                                                         &
 511                                                                                                           !faf2d,fbf2d,fab2d,fbb2d,                                                             &
 512                                                                                                           !faf2p,fbf2p,fab2p,fbb2p,                                                             &
 513                                                                                                           !faf2c,fbf2c,fab2c,fbb2c,                                                             &
 514                                                                                                            faf3,fbf3,fab3,fbb3,                                                                 &
 515                                                                                                            faf4,fbf4,fab4,fbb4,                                                                 &
 516                                                                                                            faf5,fbf5,fab5,fbb5
 517           integer,intent(in),dimension(:) :: temps_water,temps_fd,temps_crystals,                                       &
 518                                                                                  temps_snow,temps_graupel,temps_hail,                                   &
 519                                                                                  usetables
 520           real(kind=r4size),intent(in),dimension(:) :: fws_fd,fws_crystals,fws_snow,fws_graupel,fws_hail
 521
 522 ! ### Interface Vars.
 523
 524 ! ### Local Vars.
 525
 526         real(kind=r8size) :: bin_conc(number_bin)
 527         real(kind=r8size) :: ldr, kdp, cdr, ah, adp, zh
 528
 529         complex(8) :: dc_water, dc_ice, dc_wet, rhv, fa, fb, fa0, fb0, dc_wet_core
 530         complex(8) :: dc_wet_inner, dc_snow
 531         complex(8) :: sum_rhv, &
 532                    ssum_rhv
 533
 534         complex(8) :: f_a(number_bin), f_b(number_bin), &
 535                  f_a0(number_bin), f_b0(number_bin)
 536
 537         real(kind=r8size) :: a_w(7), a_column(7), a(7,number_bin)
 538         real(kind=r8size), parameter :: pi = 3.14159265D0, den_water = 1.0d0, den_ice = 0.91, den_grau0 = 0.4
 539         real(kind=r8size) :: sum_zh, sum_zv, ssum_zv, sum_ldr, sum_kdp, ssum_zh, ssum, zv, ssum_ldr, ssum_kdp,    &
 540                                                 sum_cdr, ssum_cdr, sum_ah, sum_adp, ssum_ah, ssum_adp, degree, z, x, elev, &
 541                                                 temperature, b_mass, water_mass, coef1, coef2, coef3, hail_mass, fract_mass_water, &
 542                                                 density_average, fvw, fd_mass, density_bulk, grau_mass, fract_water_crit,fract_water_scaled, &
 543                                                 fvw_core, density_core, plate_mass, dendr_mass, bulk_mass, dendr_log, density_dry, &
 544                                                 dd_dry, snow_mass, snow_log, beta, colum_mass, colum_log
 545         integer :: kb, i, itemp_w(size(temps_water)), infw_w, itemp_fd(size(temps_fd)), infw_fd(size(fws_fd)), &
 546                            itemp_g(size(temps_graupel)), infw_g(size(fws_graupel)), itemp_h(size(temps_hail)), infw_h(size(fws_hail)), &
 547                            itemp_s(size(temps_snow)), infw_s(size(fws_snow))
 548
 549 ! ### Local Vars.
 550
 551   itemp_w = 0
 552
 553
 554
 555 !**** **************************************************
 556 ! General input &
 557 !**** ********************
 558       sum_zh   =  0.0d0
 559       sum_zv   =  0.0d0
 560       sum_ldr  =  0.0d0
 561       sum_kdp  =  0.0d0
 562       ssum_zh  =  0.0d0
 563       ssum_zv  =  0.0d0
 564       ssum_ldr =  0.0d0
 565       ssum_kdp =  0.0d0
 566       sum_rhv  = (0.d0,0.d0)
 567       ssum_rhv = (0.d0,0.d0)
 568       sum_cdr = 0.0d0
 569       ssum_cdr = 0.0d0
 570       sum_ah = 0.0d0
 571       sum_adp = 0.0d0
 572       ssum_ah=0.0d0
 573       ssum_adp=0.0d0
 574
 575       degree=1.0d0/3.0d0
 576
 577 !     z = dz*(kz-1)
 578       z = zmks_1d(kz)
 579       x = dx*(kx-1)
 580
 581 ! JCS - Set elevation to 0 degrees
 582       !elev=atan(z/(x+distance))
 583           elev = 0.0d0
 584
 585 !**** ******************************************************
 586 ! Water and ice dielectric constant &
 587 !**** *************************************
 588
 589       temperature = temp-273.15d0
 590
 591       call calc_dc_water(temperature, wavelength, dc_water)
 592
 593 !*** JCS - Don't allow ice to exceed 0.0 C!
 594       call calc_dc_ice(min(temperature,0.0D0), wavelength, dc_ice)
 595
 596 !**** ******************************************************
 597 ! Water amplitude &
 598 !**** *************************
 599 ! Andrei's new change of 4.08.11                              (start)
 600
 601       call calc_orient_water(a_w)
 602
 603 ! Andrei's new change of 4.08.11                                (end)
 604
 605       do kb=1,number_bin
 606
 607          bin_conc(kb)= 0.23105d6*FF1(kb)*RORD/bin_mass(kb)
 608          b_mass      = bin_conc(kb)*bin_mass(kb)
 609
 610          if(b_mass.gt.1.0d-8) then
 611
 612            water_mass = bin_mass(kb)
 613
 614                  if(usetables(1) == 1) then
 615                itemp_w = minloc(abs(dble(temps_water)-temperature))
 616                infw_w = 1
 617                f_a(kb)  = fab1(kb,1,itemp_w(1),iwl)
 618                f_b(kb)  = fbb1(kb,1,itemp_w(1),iwl)
 619                f_a0(kb) = faf1(kb,1,itemp_w(1),iwl)
 620                f_b0(kb) = fbf1(kb,1,itemp_w(1),iwl)
 621            else
 622
 623                                       call calc_scattering_water &
 624                                                         (wavelength, water_mass,dc_water,fa,fb,fa0,fb0)
 625
 626                                                f_a(kb)  = fa
 627                                    f_b(kb)  = fb
 628                                    f_a0(kb) = fa0
 629                                              f_b0(kb) = fb0
 630             endif
 631
 632          else
 633
 634            f_a(kb)  = (0.d0,0.d0)
 635            f_b(kb)  = (0.d0,0.d0)
 636            f_a0(kb) = (0.d0,0.d0)
 637            f_b0(kb) = (0.d0,0.d0)
 638
 639          endif
 640
 641       enddo
 642 ! cycle by kb
 643
 644 ! Andrei's new change of 4.08.11                              (start)
 645       do kb=1,number_bin  ! ### (KS)
 646          do i=1,7
 647                      a(i,kb)=a_w(i)
 648                      enddo
 649       enddo
 650 ! Andrei's new change of 4.08.11                                (end)
 651
 652       call integr &
 653                         (a,bin_conc,f_a,f_b,f_a0,f_b0,zh,zv,ldr,kdp,rhv,cdr,ah,adp,ijk,kx,kz,1,number_bin)
 654
 655       coef1 = 4.0d4*(wavelength/pi)**4/ &
 656                      abs((dc_water-1.0d0)/(dc_water+2.0d0))**2
 657
 658       coef2 = 0.18d1*wavelength/pi
 659       coef3 = 8.686d-2*wavelength
 660
 661       sum_zh  = coef1*zh
 662       sum_zv  = coef1*zv
 663       sum_ldr = coef1*ldr
 664       sum_kdp = coef2*kdp
 665       sum_rhv = coef1*rhv
 666       !sum_cdr = (zh+zv-2*abs(rhv))/(zh+zv+2*abs(rhv))
 667       sum_cdr = cdr
 668       sum_ah = coef3*ah
 669       sum_adp = coef3*adp
 670
 671       call output(sum_zh,sum_zv,sum_ldr,sum_kdp,sum_rhv,sum_cdr,sum_ah,sum_adp,out1)
 672
 673 !**** ********************************************************
 674 ! Hail  amplitude &
 675 !**** *******************
 676
 677       do kb=1,number_bin
 678
 679          bin_conc(kb)=0.23105d6*FF5(kb)*RORD/bin_mass(kb)
 680
 681          b_mass=bin_conc(kb)*bin_mass(kb)
 682
 683          if(b_mass.gt.1.0d-8) then
 684
 685            hail_mass=bin_mass(kb)
 686
 687            if(FL5(kb) < 0.01d0) FL5(kb) = 0.01d0
 688
 689            fract_mass_water=FL5(kb)
 690
 691            fvw = den_ice*fract_mass_water/((1-fract_mass_water)*den_water+fract_mass_water*den_ice)
 692
 693                if(usetables(6) == 1) then
 694                  itemp_h = minloc(abs(temps_hail-temperature))
 695                  infw_h = minloc(abs(fws_hail-fvw))
 696                  f_a(kb)  = fab5(kb,infw_h(1),itemp_h(1),iwl)
 697                  f_b(kb)  = fbb5(kb,infw_h(1),itemp_h(1),iwl)
 698                  f_a0(kb) = faf5(kb,infw_h(1),itemp_h(1),iwl)
 699                  f_b0(kb) = fbf5(kb,infw_h(1),itemp_h(1),iwl)
 700                                  !if (f_a(kb)*f_b(kb)*f_a0(kb)*f_b0(kb) == 0.0d0) then
 701                  ! print *,'One of the scattering amplitudes for kb=',kb,' is 0. FIX THIS!'
 702                  !endif
 703             else
 704
 705                                      density_average=(1.0d0-fvw)*den_ice+fvw*den_water
 706   ! JCS - Although calc_dc_wet_snow uses fvw, it's calculated in the subroutine
 707   ! using den_ice, density of water, and fract_mass_water
 708                                        call calc_dc_wet_snow &
 709                                                                        (den_ice,fract_mass_water,dc_water,dc_ice,dc_wet)
 710                  call calc_scattering_hail(wavelength,hail_mass, &
 711                                                                        den_ice,fract_mass_water,dc_water,dc_ice,dc_wet,fa,fb,fa0,fb0)
 712
 713                                         f_a(kb)  = fa
 714                                         f_b(kb)  = fb
 715                                         f_a0(kb) = fa0
 716                                         f_b0(kb) = fb0
 717               endif
 718
 719 ! new change 4.08.11                                         (start)
 720                       call calc_orient(fract_mass_water,a,kb,number_bin) ! ### (KS)
 721 ! new change 4.08.11                                           (end)
 722          else
 723            f_a(kb)  = (0.d0,0.d0)
 724            f_b(kb)  = (0.d0,0.d0)
 725            f_a0(kb) = (0.d0,0.d0)
 726            f_b0(kb) = (0.d0,0.d0)
 727          endif
 728       enddo
 729 ! cycle by kb
 730
 731       call integr &
 732                   (a,bin_conc,f_a,f_b,f_a0,f_b0,zh,zv,ldr,kdp,rhv,cdr,ah,adp,ijk,kx,kz,2,number_bin)
 733
 734       coef1 = 4.0d4*(wavelength/pi)**4/ &
 735               abs((dc_water-1.0d0)/(dc_water+2.0d0))**2
 736
 737       coef2 = 0.18d1*wavelength/pi
 738
 739       ssum_zh  = coef1*zh
 740       ssum_zv  = coef1*zv
 741       ssum_ldr = coef1*ldr
 742       ssum_kdp = coef2*kdp
 743       ssum_rhv = coef1*rhv
 744       !ssum_cdr = (zh+zv-2*abs(rhv))/(zh+zv+2*abs(rhv))
 745       ssum_cdr = cdr
 746       ssum_ah = coef3*ah
 747       ssum_adp = coef3*adp
 748
 749       sum_zh  = sum_zh + ssum_zh
 750       sum_zv  = sum_zv + ssum_zv
 751       sum_ldr = sum_ldr + ssum_ldr
 752       sum_kdp = sum_kdp + ssum_kdp
 753       sum_rhv = sum_rhv + ssum_rhv
 754       sum_cdr = sum_cdr + ssum_cdr
 755       sum_ah  = sum_ah + ssum_ah
 756       sum_adp = sum_adp + ssum_adp
 757
 758       call output(ssum_zh,ssum_zv,ssum_ldr,ssum_kdp,ssum_rhv,ssum_cdr,&
 759                   ssum_ah,ssum_adp,out2)
 760
 761 !**** ********************************************************
 762 ! Freezing drops amplitude &
 763 !**** *******************
 764 ! ###################################### !
 765 ! We currently do not have FD in WRF-SBM
 766 ! ###################################### !
 767
 768
 769       do kb=1,number_bin
 770
 771          bin_conc(kb)=0.23105d6*FF1_FD(kb)*RORD/bin_mass(kb)
 772
 773          b_mass=bin_conc(kb)*bin_mass(kb)
 774
 775          if(b_mass.gt.1.0d-8) then
 776
 777            fd_mass = bin_mass(kb)
 778
 779            if(FL1_FD(kb).lt.0.01d0) FL1_FD(kb)=0.01d0
 780
 781            fract_mass_water=FL1_FD(kb)
 782
 783            density_bulk = &
 784                                 den_water*den_ice*(1.0d0-fract_mass_water)/ &
 785                                 (den_water*(1.0d0-fract_mass_water)+ &
 786                                 den_ice*fract_mass_water)
 787
 788            fvw = den_ice*fract_mass_water/((1-fract_mass_water)*den_water+fract_mass_water*den_ice)
 789
 790                         if(usetables(2) == 1) then
 791                             !itemp_fd = minloc(abs(temps_fd-temperature))
 792                             !infw_fd = minloc(abs(fws_fd-fvw))
 793                             !f_a(kb)  = fab1fd(kb,infw_fd(1),itemp_fd(1),iwl)
 794                             !f_b(kb)  = fbb1fd(kb,infw_fd(1),itemp_fd(1),iwl)
 795                             !f_a0(kb) = faf1fd(kb,infw_fd(1),itemp_fd(1),iwl)
 796                             !f_b0(kb) = fbf1fd(kb,infw_fd(1),itemp_fd(1),iwl)
 797                         else
 798
 799                                 density_average=(1.0d0-fvw)*den_ice+fvw*den_water
 800
 801                                 call calc_dc_wet_snow &
 802                                                         (den_ice,fract_mass_water,dc_water,dc_ice,dc_wet)
 803
 804 ! JCS -- Needed to modify argument list to pass in dc_wet for tmatrix
 805 ! calculations
 806                                 call calc_scattering_fd &
 807                                                 (wavelength,fd_mass, &
 808                                                 density_average,fract_mass_water,dc_water,dc_ice,dc_wet, &
 809                                                 fa,fb,fa0,fb0)
 810
 811                                 f_a(kb)  = fa
 812                                 f_b(kb)  = fb
 813                                 f_a0(kb) = fa0
 814                                 f_b0(kb) = fb0
 815                         endif
 816
 817                         call calc_orient(fract_mass_water,a,kb,number_bin) ! ### (KS)
 818
 819          else
 820
 821            f_a(kb)  = (0.d0,0.d0)
 822            f_b(kb)  = (0.d0,0.d0)
 823            f_a0(kb) = (0.d0,0.d0)
 824            f_b0(kb) = (0.d0,0.d0)
 825
 826          endif
 827
 828       enddo
 829 ! cycle by kb
 830
 831       call integr &
 832                 (a,bin_conc,f_a,f_b,f_a0,f_b0,zh,zv,ldr,kdp,rhv,cdr,ah,adp,ijk,kx,kz,3,number_bin)
 833
 834       coef1 = 4.0d4*(wavelength/pi)**4/ &
 835               abs((dc_water-1.0d0)/(dc_water+2.0d0))**2
 836
 837       coef2 = 0.18d1*wavelength/pi
 838
 839       ssum_zh  = coef1*zh
 840       ssum_zv  = coef1*zv
 841       ssum_ldr = coef1*ldr
 842       ssum_kdp = coef2*kdp
 843       ssum_rhv = coef1*rhv
 844       !ssum_cdr = (zh+zv-2*abs(rhv))/(zh+zv+2*abs(rhv))
 845       ssum_cdr = cdr
 846       ssum_ah  = coef3*ah
 847       ssum_adp = coef3*adp
 848
 849       sum_zh  = sum_zh +ssum_zh
 850       sum_zv  = sum_zv +ssum_zv
 851       sum_ldr = sum_ldr+ssum_ldr
 852       sum_kdp = sum_kdp+ssum_kdp
 853       sum_rhv = sum_rhv+ssum_rhv
 854       sum_ah  = sum_ah +ssum_ah
 855       sum_adp = sum_adp+ssum_adp
 856       sum_cdr = sum_cdr+ssum_cdr
 857
 858       call output(ssum_zh,ssum_zv,ssum_ldr,ssum_kdp,ssum_rhv,ssum_cdr,&
 859                                 ssum_ah,ssum_adp,out3)
 860
 861
 862 !**** ***************************************************************
 863 ! Graupel  amplitude &
 864 !**** ********************
 865
 866       do kb=1,number_bin
 867
 868          bin_conc(kb)=0.23105d6*FF4(kb)*RORD/bin_mass(kb)
 869
 870          b_mass=bin_conc(kb)*bin_mass(kb)
 871
 872          if(b_mass.gt.1.d-8) then
 873
 874            grau_mass= bin_mass(kb)
 875            fract_mass_water=FL4(kb)
 876
 877 ! new change 4.08.11                                         (start)
 878            call calc_orient(fract_mass_water,a,kb,number_bin) ! ### (KS)
 879 ! new change 4.08.11                                           (end)
 880 ! JCS - I don't know where the following equation comes from. Regardless,
 881 ! we'll use it to define the maximum fractional water that a particle can have
 882 ! before it begins to have a water coating.
 883            fract_water_crit = den_water*(den_ice-den_grau0)/ &
 884                                                         (den_water*(den_ice-den_grau0)+den_ice*den_grau0)
 885
 886            density_bulk =  &
 887                                 den_water*den_grau0*(1.0d0-fract_mass_water)/ &
 888                                 (den_water*(1.0d0-fract_mass_water)+ &
 889                                 den_grau0*fract_mass_water)
 890
 891 ! JCS - Calculate the average density of the graupel varying between
 892 ! den_grau0 and den_water based upon fractional volume of water. The average
 893 ! density needs the fractional volume of water (not fractional mass of water)!
 894            fvw = den_grau0*fract_mass_water/((1-fract_mass_water)*den_water+fract_mass_water*den_grau0)
 895
 896                    if(usetables(5) == 1) then
 897                itemp_g = minloc(abs(temps_graupel-temperature))
 898                infw_g = minloc(abs(fws_graupel-fvw))
 899                f_a(kb)  = fab4(kb,infw_g(1),itemp_g(1),iwl)
 900                f_b(kb)  = fbb4(kb,infw_g(1),itemp_g(1),iwl)
 901                f_a0(kb) = faf4(kb,infw_g(1),itemp_g(1),iwl)
 902                f_b0(kb) = fbf4(kb,infw_g(1),itemp_g(1),iwl)
 903            else
 904
 905                                 density_average=(1.0d0-fvw)*den_grau0+fvw*den_water
 906 ! JCS - DC_wet will be between the DC of a dry graupel particle (with density of
 907 ! den_grau0) and the DC of a water drop).
 908                         call calc_dc_wet_snow &
 909                                                                         (den_grau0,fract_mass_water, &
 910                                                                         dc_water,dc_ice,dc_wet)
 911
 912                         if(fract_mass_water.lt.fract_water_crit) then
 913 ! JCS - Model graupel as spongy and use the dc_wet obtained above
 914                         call calc_scattering_grau1 &
 915                                                                         (wavelength,grau_mass,density_average, &
 916                                                                         fract_mass_water,dc_wet,fa,fb,fa0,fb0)
 917                         else
 918 ! JCS - fract_water_scaled is the fractional water in the soaked ice core. The
 919 ! rest of the water is going to coat the particle. We need to find the DC of the
 920 ! soaked inner core, which means we'll need the fvw of only the core (ignore the
 921 ! excess water that'll coat the particle).
 922 ! Barry/Kobby Correction  FRACT_CRIT_WATER
 923                         fract_water_scaled = fract_water_crit/(1-fract_mass_water+fract_water_crit)
 924                         fvw_core = den_grau0*fract_water_scaled/((1-fract_water_scaled)*den_water+ &
 925                                fract_water_scaled*den_grau0)
 926 ! JCS - density_core is the density of the interior of the soaked particle.
 927 ! We'll use fvw_core to find what's essentially the critical density
 928                         density_core=(1.d0-fvw_core)*den_grau0+fvw_core*den_water
 929 ! JCS - calculate dc_wet_inner as the dielectric constant of the soaked inner
 930 ! spongy core. This will only be used for two-layer calculations. The
 931 ! dc_wet, which is for the entire particle, is used for T-matrix calculations
 932 ! since the T-matrix code assumes a homogeneous mixture at this time.
 933                         call calc_dc_wet_snow &
 934                                                                 (den_grau0,fract_water_scaled, &
 935                                                                 dc_water,dc_ice,dc_wet_core)
 936 ! Model using two-layer if the particle is Rayleigh-sized; if it's larger, then
 937 ! it'll be modeled as a homogeneous mixture using the t_matrix.F90 subroutine.
 938                                         call calc_scattering_grau2 &
 939                                                                                 (wavelength,grau_mass, &
 940                                                                                 density_core,fract_mass_water,fract_water_crit,dc_water,dc_ice,dc_wet_core, &
 941                                                                                 dc_wet,fa,fb,fa0,fb0)
 942                         endif
 943                                 ! in case fract_mass_water.lt.fract_water_crit
 944
 945            f_a(kb)  = fa
 946            f_b(kb)  = fb
 947            f_a0(kb) = fa0
 948            f_b0(kb) = fb0
 949
 950                 endif
 951
 952
 953          else
 954         ! in case b_mass.le.1.d-8
 955
 956                 f_a(kb)  = (0.d0,0.d0)
 957                 f_b(kb)  = (0.d0,0.d0)
 958                 f_a0(kb) = (0.d0,0.d0)
 959                 f_b0(kb) = (0.d0,0.d0)
 960
 961          endif
 962   enddo
 963 ! cycle by kb
 964
 965       call integr &
 966                         (a,bin_conc,f_a,f_b,f_a0,f_b0,zh,zv,ldr,kdp,rhv,cdr,ah,adp,ijk,kx,kz,4,number_bin)
 967
 968       coef1 = 4.0d4*(wavelength/pi)**4/ &
 969               abs((dc_water-1.0d0)/(dc_water+2.0d0))**2
 970
 971       coef2 = 0.18d1*wavelength/pi
 972
 973       ssum_zh  = coef1*zh
 974       ssum_zv  = coef1*zv
 975       ssum_ldr = coef1*ldr
 976       ssum_kdp = coef2*kdp
 977       ssum_rhv = coef1*rhv
 978       ssum_cdr = cdr
 979       ssum_ah = coef3*ah
 980       ssum_adp = coef3*adp
 981
 982       sum_zh  = sum_zh +ssum_zh
 983       sum_zv  = sum_zv +ssum_zv
 984       sum_ldr = sum_ldr+ssum_ldr
 985       sum_kdp = sum_kdp+ssum_kdp
 986       sum_rhv = sum_rhv+ssum_rhv
 987       sum_cdr = sum_cdr+ssum_cdr
 988       sum_ah = sum_ah+ssum_ah
 989       sum_adp = sum_adp+ssum_adp
 990
 991       call output(ssum_zh,ssum_zv,ssum_ldr,ssum_kdp,ssum_rhv,ssum_cdr,&
 992                                 ssum_ah,ssum_adp,out4)
 993
 994 !**** ********************************************************
 995 ! Plate  amplitude   * &
 996 !**** ******************
 997
 998       do kb=1,number_bin
 999
1000          bin_conc(kb)=0.23105d6*FF2(kb,2)*RORD/bin_mass(kb)
1001
1002          b_mass=bin_conc(kb)*bin_mass(kb)
1003
1004          if(b_mass.gt.1.d-8) then
1005
1006            plate_mass= bin_mass(kb)
1007            fract_mass_water=0.0d0
1008
1009                    ! new change 4.08.11                                         (start)
1010            call calc_orient(fract_mass_water,a,kb,number_bin) ! ### (KS)
1011                         ! new change 4.08.11                                           (end)
1012
1013            call calc_rayleigh_plate &
1014                                                         (wavelength,plate_mass, &
1015                                                         den_ice,fract_mass_water,dc_water,dc_ice, &
1016                                                         fa,fb,fa0,fb0)
1017
1018            f_a(kb)  = fa
1019            f_b(kb)  = fb
1020            f_a0(kb) = fa0
1021            f_b0(kb) = fb0
1022
1023          else
1024
1025            f_a(kb)  = (0.d0,0.d0)
1026            f_b(kb)  = (0.d0,0.d0)
1027            f_a0(kb) = (0.d0,0.d0)
1028            f_b0(kb) = (0.d0,0.d0)
1029
1030          endif
1031
1032       enddo
1033
1034       call integr &
1035                                 (a,bin_conc,f_a,f_b,f_a0,f_b0,zh,zv,ldr,kdp,rhv,cdr,ah,adp,ijk,kx,kz,5,number_bin)
1036
1037       coef1 = 4.0d4*(wavelength/pi)**4/ &
1038               abs((dc_water-1.0d0)/(dc_water+2.0d0))**2
1039
1040       coef2 = 0.18d1*wavelength/pi
1041
1042       ssum_zh  = coef1*zh
1043       ssum_zv  = coef1*zv
1044       ssum_ldr = coef1*ldr
1045       ssum_kdp = coef2*kdp
1046       ssum_rhv = coef1*rhv
1047       ssum_cdr = cdr
1048       ssum_ah = coef3*ah
1049       ssum_adp = coef3*adp
1050
1051       sum_zh  = sum_zh +ssum_zh
1052       sum_zv  = sum_zv +ssum_zv
1053       sum_ldr = sum_ldr+ssum_ldr
1054       sum_kdp = sum_kdp+ssum_kdp
1055       sum_rhv = sum_rhv+ssum_rhv
1056       sum_cdr = sum_cdr+ssum_cdr
1057       sum_ah = sum_ah+ssum_ah
1058       sum_adp = sum_adp+ssum_adp
1059
1060       call output(ssum_zh,ssum_zv,ssum_ldr,ssum_kdp,ssum_rhv,ssum_cdr,&
1061                                 ssum_ah,ssum_adp,out5)
1062
1063 !**** ********************************************************
1064 ! Dendrit amplitude   * &
1065 !**** ********************
1066
1067       do kb=1,number_bin
1068
1069          bin_conc(kb)=0.23105d6*FF2(kb,3)*RORD/bin_mass(kb)
1070
1071          b_mass=bin_conc(kb)*bin_mass(kb)
1072
1073          if(b_mass.gt.1.d-8) then
1074
1075            dendr_mass= bin_mass(kb)
1076            fract_mass_water=0.0d0
1077
1078 ! new change 4.08.11                                         (start)
1079            call calc_orient(fract_mass_water,a,kb,number_bin) ! ### (KS)
1080 ! new change 4.08.11                                           (end)
1081
1082            bulk_mass = (1.0d0-fract_mass_water)*dendr_mass
1083            dendr_log =  log10(bulk_mass)
1084
1085            call INTERPOL &
1086                                 (number_bin,bin_log,tab_dendr,dendr_log,density_bulk)
1087
1088            dendr_log  = log10(dendr_mass)
1089
1090            call INTERPOL &
1091                                 (number_bin,bin_log,tab_dendr,dendr_log,density_dry)
1092
1093            dd_dry = 1.d1*(dendr_mass/density_dry)**degree
1094
1095            call calc_dc_wet_snow &
1096                                 (density_bulk,fract_mass_water,dc_water,dc_ice,dc_wet)
1097
1098            call calc_rayleigh_dendr &
1099                                 (wavelength,dendr_mass, &
1100                                 density_bulk,fract_mass_water,dd_dry,dc_wet, &
1101                                 fa,fb,fa0,fb0,ijk,kx,kz,kb)
1102
1103            f_a(kb)  = fa
1104            f_b(kb)  = fb
1105            f_a0(kb) = fa0
1106            f_b0(kb) = fb0
1107
1108 ! in case b_mass.gt.1.d-8
1109
1110          else
1111
1112 ! in case b_mass.le.1.d-8
1113
1114            f_a(kb)  = (0.d0,0.d0)
1115            f_b(kb)  = (0.d0,0.d0)
1116            f_a0(kb) = (0.d0,0.d0)
1117            f_b0(kb) = (0.d0,0.d0)
1118
1119          endif
1120
1121       enddo
1122
1123       call integr &
1124                         (a,bin_conc,f_a,f_b,f_a0,f_b0,zh,zv,ldr,kdp,rhv,cdr,ah,adp,ijk,kx,kz,6,number_bin)
1125
1126       coef1 = 4.0d4*(wavelength/pi)**4/ &
1127               abs((dc_water-1.0d0)/(dc_water+2.0d0))**2
1128
1129       coef2 = 0.18d1*wavelength/pi
1130
1131       ssum_zh  = coef1*zh
1132       ssum_zv  = coef1*zv
1133       ssum_ldr = coef1*ldr
1134       ssum_kdp = coef2*kdp
1135       ssum_rhv = coef1*rhv
1136       ssum_cdr = cdr
1137       ssum_ah = coef3*ah
1138       ssum_adp = coef3*adp
1139
1140       sum_zh  = sum_zh +ssum_zh
1141       sum_zv  = sum_zv +ssum_zv
1142       sum_ldr = sum_ldr+ssum_ldr
1143       sum_kdp = sum_kdp+ssum_kdp
1144       sum_rhv = sum_rhv+ssum_rhv
1145       sum_cdr = sum_cdr+ssum_cdr
1146       sum_ah = sum_ah+ssum_ah
1147       sum_adp = sum_adp+ssum_adp
1148
1149       call output(ssum_zh,ssum_zv,ssum_ldr,ssum_kdp,ssum_rhv,ssum_cdr,&
1150                                   ssum_ah,ssum_adp,out6)
1151
1152
1153 !**** ********************************************************
1154 ! Snow flakes  amplitude   * &
1155 !**** ***********************
1156
1157       do kb=1,number_bin
1158
1159          bin_conc(kb)=0.23105d6*FF3(kb)*RORD/bin_mass(kb)
1160
1161          b_mass=bin_conc(kb)*bin_mass(kb)
1162
1163          if(b_mass.gt.1.d-8) then
1164
1165            snow_mass        = bin_mass(kb)
1166            fract_mass_water = FL3(kb)
1167
1168 ! new change 4.08.11                                         (start)
1169            call calc_orient(fract_mass_water,a,kb,number_bin) ! ### (KS)
1170 ! new change 4.08.11                                           (end)
1171
1172            density_bulk = bulk(kb)
1173
1174            snow_log  = log10(snow_mass)
1175
1176            call INTERPOL &
1177                                 (number_bin,bin_log,tab_snow,snow_log,density_dry)
1178
1179                    fvw=density_dry*fract_mass_water/((1-fract_mass_water)*den_water+fract_mass_water*density_dry)
1180
1181            if(usetables(4) == 1) then
1182                itemp_s = minloc(abs(temps_snow-temperature))
1183                infw_s = minloc(abs(fws_snow-fvw))
1184                f_a(kb)  = fab3(kb,infw_s(1),itemp_s(1),iwl)
1185                f_b(kb)  = fbb3(kb,infw_s(1),itemp_s(1),iwl)
1186                f_a0(kb) = faf3(kb,infw_s(1),itemp_s(1),iwl)
1187                f_b0(kb) = fbf3(kb,infw_s(1),itemp_s(1),iwl)
1188                !if (f_a(kb)*f_b(kb)*f_a0(kb)*f_b0(kb) == 0.0d0) then
1189                !  print *,'One of the scattering amplitudes (SNOW) for kb=',kb,' is 0. FIX THIS!'
1190                !endif
1191            else
1192
1193                                 density_average=(1.0d0-fvw)*density_dry+fvw*den_water
1194 ! JCS - Although calc_dc_wet_snow uses fvw, it's calculated in the subroutine
1195 ! using den_ice, density of water, and fract_mass_water
1196                                 call calc_dc_dry_snow(density_dry,dc_ice,dc_snow)
1197                                 call calc_dc_wet_snow(density_average,fract_mass_water,dc_water,dc_ice,dc_wet)
1198                                 call calc_scattering_snow(wavelength,snow_mass,density_dry,density_average, &
1199                                           fract_mass_water,dc_water,dc_snow,dc_wet,fa,fb,fa0,fb0)
1200                                 f_a(kb)  = fa
1201                                 f_b(kb)  = fb
1202                                 f_a0(kb) = fa0
1203                                 f_b0(kb) = fb0
1204            endif
1205 ! new change 4.08.11                                         (start)
1206            call calc_orient(fract_mass_water,a,kb,number_bin)  ! ### [KS]
1207 ! new change 4.08.11                                           (end)
1208          else
1209
1210            f_a(kb)  = (0.d0,0.d0)
1211            f_b(kb)  = (0.d0,0.d0)
1212            f_a0(kb) = (0.d0,0.d0)
1213            f_b0(kb) = (0.d0,0.d0)
1214
1215          endif
1216       enddo
1217 ! cycle by kb
1218
1219       call integr &
1220                         (a,bin_conc,f_a,f_b,f_a0,f_b0,zh,zv,ldr,kdp,rhv,cdr,ah,adp,ijk,kx,kz,7,number_bin)
1221
1222       coef1 = 4.0d4*(wavelength/pi)**4/ &
1223               abs((dc_water-1.0d0)/(dc_water+2.0d0))**2
1224
1225       coef2 = 0.18d1*wavelength/pi
1226
1227       ssum_zh  = coef1*zh
1228       ssum_zv  = coef1*zv
1229       ssum_ldr = coef1*ldr
1230       ssum_kdp = coef2*kdp
1231       ssum_rhv = coef1*rhv
1232       ssum_cdr = cdr
1233       ssum_ah = coef3*ah
1234       ssum_adp = coef3*adp
1235
1236       sum_zh  = sum_zh +ssum_zh
1237       sum_zv  = sum_zv +ssum_zv
1238       sum_ldr = sum_ldr+ssum_ldr
1239       sum_kdp = sum_kdp+ssum_kdp
1240       sum_rhv = sum_rhv+ssum_rhv
1241       sum_cdr = sum_cdr+ssum_cdr
1242       sum_ah = sum_ah+ssum_ah
1243       sum_adp = sum_adp+ssum_adp
1244
1245       call output(ssum_zh,ssum_zv,ssum_ldr,ssum_kdp,ssum_rhv,ssum_cdr,&
1246                           ssum_ah,ssum_adp,out7)
1247
1248 !**** ********************************************************
1249 ! Column  amplitude   * &
1250 !**** *******************
1251 ! Andrei's new change of 4.08.11                             (start)
1252
1253       fract_mass_water = 0.0d0
1254
1255       beta = elev
1256
1257       call calc_orient_colum(beta,a_column)
1258
1259 ! Andrei's new change of 4.08.11                                (end)
1260
1261       do kb=1,number_bin
1262
1263          bin_conc(kb)=0.23105d6*FF2(kb,1)*RORD/bin_mass(kb)
1264
1265          b_mass=bin_conc(kb)*bin_mass(kb)
1266
1267          if(b_mass.gt.1.d-8) then
1268
1269            colum_mass= bin_mass(kb)
1270            bulk_mass = (1.0d0-fract_mass_water)*colum_mass
1271            colum_log = log10(bulk_mass)
1272
1273            call INTERPOL &
1274                                 (number_bin,bin_log,tab_colum,colum_log,density_bulk)
1275
1276            colum_log = log10(colum_mass)
1277
1278            call INTERPOL &
1279                                 (number_bin,bin_log,tab_colum,colum_log,density_dry)
1280
1281            dd_dry = 1.d1*(colum_mass/density_dry)**degree
1282
1283            call calc_dc_wet_snow &
1284                                 (density_bulk,fract_mass_water,dc_water,dc_ice,dc_wet)
1285
1286            call calc_rayleigh_colum &
1287                                 (wavelength,colum_mass,density_bulk,fract_mass_water, &
1288                                 dd_dry,dc_water,dc_wet,fa,fb,fa0,fb0)
1289
1290            f_a(kb)  = fa
1291            f_b(kb)  = fb
1292            f_a0(kb) = fa0
1293            f_b0(kb) = fb0
1294
1295          else
1296
1297            f_a(kb)  = (0.d0,0.d0)
1298            f_b(kb)  = (0.d0,0.d0)
1299            f_a0(kb) = (0.d0,0.d0)
1300            f_b0(kb) = (0.d0,0.d0)
1301
1302          endif
1303       enddo
1304 ! cycle by kb
1305
1306 ! Andrei's new change of 4.08.11                              (start)
1307
1308       do kb=1,number_bin ! ### (KS)
1309          do i=1,7
1310             a(i,kb)=a_column(i)
1311          enddo
1312       enddo
1313
1314 ! Andrei's new change of 4.08.11                                (end)
1315
1316       call integr &
1317                 (a,bin_conc,f_a,f_b,f_a0,f_b0,zh,zv,ldr,kdp,rhv,cdr,ah,adp,ijk,kx,kz,8,number_bin)
1318
1319       coef1 = 4.0d4*(wavelength/pi)**4/ &
1320                abs((dc_water-1.0d0)/(dc_water+2.0d0))**2
1321
1322       coef2 = 0.18d1*wavelength/pi
1323
1324       ssum_zh  = coef1*zh
1325       ssum_zv  = coef1*zv
1326       ssum_ldr = coef1*ldr
1327       ssum_kdp = coef2*kdp
1328       ssum_rhv = coef1*rhv
1329       ssum_cdr = cdr
1330       ssum_ah = coef3*ah
1331       ssum_adp = coef3*adp
1332
1333       sum_zh  = sum_zh +ssum_zh
1334       sum_zv  = sum_zv +ssum_zv
1335       sum_ldr = sum_ldr+ssum_ldr
1336       sum_kdp = sum_kdp+ssum_kdp
1337       sum_rhv = sum_rhv+ssum_rhv
1338       sum_cdr = sum_cdr+ssum_cdr
1339       sum_ah = sum_ah+ssum_ah
1340       sum_adp = sum_adp+ssum_adp
1341
1342       call output(ssum_zh,ssum_zv,ssum_ldr,ssum_kdp,ssum_rhv,ssum_cdr,&
1343                   ssum_ah,ssum_adp,out8)
1344
1345       call output(sum_zh,sum_zv,sum_ldr,sum_kdp,sum_rhv,sum_cdr,sum_ah,&
1346                   sum_adp,out9)
1347
1348       return
1349       end subroutine polar_hucm
1350
1351 ! subroutine polar_hucm &
1352
1353 !**** ************************************************************** &
1354 !**** **************************************************************
1355
1356       SUBROUTINE INTERPOL(NH, H_TAB, X_TAB, H, X)
1357
1358       !IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1359            implicit none
1360 ! ### Interface
1361           integer :: NH
1362       double precision :: H_TAB(NH), X_TAB(NH)
1363           double precision :: H, X
1364 ! ### Interface
1365           integer :: I, J
1366
1367       IF(H.LT.H_TAB(1)) THEN
1368
1369          X=X_TAB(1)
1370
1371          RETURN
1372
1373       ENDIF
1374
1375       IF(H.GT.H_TAB(NH)) THEN
1376
1377          X=X_TAB(NH)
1378
1379          RETURN
1380
1381       ENDIF
1382
1383       DO I=2,NH
1384
1385          IF(H.LT.H_TAB(I)) THEN
1386
1387             J=I-1
1388             X=X_TAB(J)+(X_TAB(I)-X_TAB(J))/ &
1389            (H_TAB(I)-H_TAB(J))*(H-H_TAB(J))
1390
1391             RETURN
1392
1393          ENDIF
1394
1395       ENDDO
1396
1397       RETURN
1398       END SUBROUTINE INTERPOL
1399
1400 ! SUBROUTINE INTERPOL &
1401
1402 !**** ******************************************************************** &
1403 !**** ********************************************************************
1404
1405       subroutine calc_dc_water(temp,wl,dc_water)
1406
1407       !implicit double precision (a-h,o-z)
1408       implicit none
1409 ! ### Interface
1410           complex(8) :: dc_water
1411           double precision :: temp,wl
1412 ! ### Interface
1413           double precision :: tt_rad, eps0, epsinf, alfa, wl0, rat, si, co, dc_re, dc_im
1414       double precision, parameter :: pi = 3.14159265D0, sig = 12.5664d8
1415
1416       tt_rad     = temp-25.0d0
1417       eps0   = 78.54d0*(1.d0-4.579d-3*tt_rad+1.19d-5*tt_rad**2-2.8d-8*tt_rad**3)
1418       epsinf = 5.27137d0+0.021647*temp-0.00131198*temp**2
1419       tt_rad     = temp+273.0d0
1420       alfa   = -16.8129d0/tt_rad+0.0609265d0
1421       wl0    = 0.00033836d0*exp(2513.98d0/tt_rad)
1422       rat    = (wl0/wl)**(1.0d0-alfa)
1423       si     = sin(0.5d0*alfa*pi)
1424       co     = cos(0.5d0*alfa*pi)
1425
1426       dc_re  = epsinf+(eps0-epsinf)*(1.0d0+rat*si)/ &
1427                       (1.0d0+2.0d0*rat*si+rat**2)
1428       dc_im  = (eps0-epsinf)*rat*co/(1.0d0+2.0d0*rat*si+rat**2)+ &
1429                sig*wl/18.8496d10
1430
1431       dc_water = cmplx(dc_re,dc_im)
1432
1433       return
1434       end subroutine calc_dc_water
1435
1436 ! subroutine calc_dc_water &
1437
1438 !**** ****************************************************************** &
1439 !**** ******************************************************************
1440
1441       subroutine calc_dc_ice(temp, wl, dc_ice)
1442
1443       !implicit double precision (a-h,o-z)
1444       implicit none
1445 ! ### Interface
1446           complex(8) :: dc_ice
1447           double precision :: temp, wl
1448 ! ### Interface
1449           double precision :: eps0, alfa, tt_rad, wl0, sig, rat, si, co, dc_re, dc_im
1450       double precision,parameter :: pi = 3.14159265D0, epsinf = 3.168d0
1451
1452       eps0   = 203.168d0+2.5d0*temp+0.15d0*temp**2
1453       alfa   = 0.288d0+0.0052d0*temp+0.00023d0*temp**2
1454       tt_rad     = temp+273.0d0
1455       wl0    = 0.0009990288d0*exp(6.6435d3/tt_rad)
1456       sig    = 1.26d0*exp(-6.2912d3/tt_rad)
1457       rat    = (wl0/wl)**(1.0d0-alfa)
1458       si     = sin(0.5d0*alfa*pi)
1459       co     = cos(0.5d0*alfa*pi)
1460
1461       dc_re  = epsinf+(eps0-epsinf)*(1.0d0+rat*si)/ &
1462                       (1.0d0+2.0d0*rat*si+rat**2)
1463       dc_im  = (eps0-epsinf)*rat*co/(1.0d0+2.0d0*rat*si+rat**2)+ &
1464                sig*wl/18.8496d10
1465
1466       dc_ice = cmplx(dc_re,dc_im)
1467
1468       return
1469       end  subroutine calc_dc_ice
1470
1471 ! subroutine calc_dc_ice &
1472
1473 !**** ****************************************************************** &
1474 !**** ******************************************************************
1475
1476       subroutine calc_dc_dry_snow(den_bulk,dc_ice,dc_snow)
1477
1478       !implicit double precision (a-h,o-z)
1479        implicit none
1480 ! ### INterface
1481           double precision :: den_bulk
1482       complex(8) :: dc_ice, dc_snow, ratc ! ### [KS] : complex(8)
1483 ! ### Interface
1484
1485       double precision,parameter :: den_ice = 0.91d0
1486           double precision :: rat
1487
1488       rat     =  den_bulk/den_ice
1489       ratc    = (dc_ice-1.0d0)/(dc_ice+2.0d0)
1490
1491       dc_snow =  1.0d0+3.0d0*rat*ratc/(1.0d0-rat*ratc)
1492
1493       return
1494       end subroutine calc_dc_dry_snow
1495
1496 ! subroutine calc_dc_dry_snow &
1497
1498 !**** ****************************************************************** &
1499 !**** ******************************************************************
1500
1501       subroutine calc_dc_wet_snow &
1502                                         (den_bulk,fract_mass_water,dc_water,dc_ice,dc_wet_snow)
1503
1504       !implicit double precision (a-h,o-z)
1505           implicit none
1506 ! ### Interface
1507       complex(8) :: dc_water,dc_ice,dc_dry_snow,dc_wet_snow1, &  ! ### [KS] : complex(8)
1508                      dc_wet_snow2,ratc,dc_wet_snow
1509           double precision :: den_bulk, fract_mass_water
1510 ! ### Interface
1511
1512       double precision,parameter :: den_water = 1.0d0, den_ice = 0.91d0
1513           double precision :: rat, fract_volume_water, t
1514
1515
1516       rat         =  den_bulk/den_ice
1517       ratc        = (dc_ice-1.0d0)/(dc_ice+2.0d0)
1518       dc_dry_snow =  1.0d0+3.0d0*rat*ratc/(1.0d0-rat*ratc)
1519 ! den_bulk should be the density of DRY snow, graupel, hail, etc.!
1520       rat                = den_bulk/den_water
1521       fract_volume_water = 1.0d0-(1.0d0-fract_mass_water)/ &
1522                           (1.0d0+fract_mass_water*(rat-1.0d0))
1523
1524       ratc = (dc_water-dc_dry_snow)/(dc_water+2.0d0*dc_dry_snow)
1525
1526       dc_wet_snow1 = dc_dry_snow* &
1527                     (1.0d0+3.0d0*fract_volume_water*ratc/ &
1528                     (1.0d0-fract_volume_water*ratc))
1529
1530       ratc = (dc_dry_snow-dc_water)/(dc_dry_snow+2.0d0*dc_water)
1531
1532       dc_wet_snow2 = dc_water* &
1533                     (1.0d0+3.0d0*(1.0d0-fract_volume_water)*ratc/ &
1534                     (1.0d0-(1.0d0-fract_volume_water)*ratc))
1535 ! new change 18.01.09                                         (start)
1536       if(fract_volume_water.gt.1.0d-10) then
1537          t=erf((1.0d0-fract_volume_water)/fract_volume_water-0.2d0)
1538       else
1539          t=1.0d0
1540       endif
1541 ! new change 18.01.09                                           (end)
1542
1543       dc_wet_snow = 0.5d0*((1.0d0+t)*dc_wet_snow1+ &
1544                            (1.0d0-t)*dc_wet_snow2)
1545
1546       return
1547       end subroutine calc_dc_wet_snow
1548
1549 ! subroutine calc_dc_wet_snow &
1550
1551 !**** *************************************************************** &
1552 !**** ***************************************************************
1553
1554       subroutine calc_scattering_water(wl,water_mass,dc,f_a,f_b,f_a0,f_b0)
1555
1556       !implicit double precision (a-h,o-z)
1557            implicit none
1558
1559            intrinsic DCONJG
1560 ! ### Interface
1561            double precision :: wl, water_mass
1562            complex(8) :: dc, f_a, f_b, f_a0, f_b0
1563 ! ### Interface
1564
1565       double precision,parameter :: pi = 3.14159265D0, den_water = 1.0d0
1566           double precision :: degree, dd, aspect, rp, angle, ff, ff2, shape_a, shape_b, tmp
1567
1568       degree=1.0d0/3.0d0
1569
1570       dd  = 1.d1*((6.0D0/pi)*water_mass/den_water)**degree
1571
1572       if(dd.lt.1.0d1) then
1573         aspect = 0.9951d0+0.0251d0*dd-0.03644*dd**2+ &
1574                  0.005303*dd**3-0.0002492*dd**4
1575       else
1576         aspect = 0.4131d0
1577       endif
1578
1579       if(aspect.eq.1.0d0) aspect=0.9999d0
1580
1581       rp = dd*abs(sqrt(dc))/(1.d1*wl)
1582
1583       if(rp.gt.0.1d0) then
1584
1585          angle = 1.8d2
1586
1587 !         call t_matrix(dd,wl,dc,aspect,angle,f_a,f_b,6,'water') ! [KS] >> This is not linked as we use lookup tables
1588
1589          angle = 0.0d0
1590
1591 !        call t_matrix(dd,wl,dc,aspect,angle,f_a0,f_b0,6,'water') ! [KS] >> This is not linked as we use lookup tables
1592          f_b0 = -DCONJG(f_b0)
1593          f_a0 = DCONJG(f_a0)
1594
1595       else
1596
1597          ff      = sqrt((1.0d0/aspect)**2-1.0d0)
1598          ff2     = ff**2
1599          shape_a = ((1.0+ff2)/ff2)*(1.0d0-atan(ff)/ff)
1600          shape_b = 0.5d0*(1.0d0-shape_a)
1601          tmp     = pi**2*dd**3/(6.0d2*wl**2)
1602          f_a0     = tmp/(shape_a+1.0d0/(dc-1.0d0))
1603          f_b0     = tmp/(shape_b+1.0d0/(dc-1.0d0))
1604          f_a    = dconjg(f_a0)
1605          f_b    = dconjg(f_b0)
1606
1607       endif
1608
1609       return
1610       end subroutine calc_scattering_water
1611
1612 ! subroutine calc_scattering_water &
1613
1614 !**** ******************************************************************** &
1615 !**** ********************************************************************
1616
1617       subroutine calc_scattering_hail &
1618                         (wl, hail_mass, den_bulk,fract_mass_water,dc_water,dc_hail,dc_wet, &
1619                          f_a,f_b,f_a0,f_b0)
1620
1621           USE scatt_tables,ONLY:twolayer_hail,rpquada,usequad
1622       !USE t_matrix2_quad_mod             ! ### [KS] : this is not linked since we use look-up-tables
1623       !USE t_matrix2_double_mod           ! ### [KS] : this is not linked since we use look-up-tables
1624
1625           implicit none
1626
1627           intrinsic DCONJG
1628      ! implicit double precision (a-h,o-z)
1629       double precision :: degree, factor, hail_mass, fvw, rpquad
1630       double precision :: wl, den_bulk, fract_mass_water, aspect_melt, aspect_dry, aspect, aspect2
1631       double precision :: dd, dd_dmelt, dd_dry, dd2,dd1, dd_melt, dcore, rp, angle
1632       double precision :: ff, ff2, shape2_a, shape2_b, shape1_a, shape1_b, psi, tmp
1633
1634       complex(8) :: dc_water,dc_hail,dc_wet,num,denum,f_a,f_b,f_a0,f_b0
1635
1636       double precision, parameter :: pi = 3.14159265D0, den_water = 1.0d0, den0 = 0.91d0
1637
1638       degree=1.0d0/3.0d0
1639
1640       dd_dry = 1.d1*((6.0D0/pi)*hail_mass/den0)**degree
1641
1642       if(dd_dry.lt.1.0d1) then
1643          aspect_dry = 1.0d0-2.0d-2*dd_dry
1644       else
1645          aspect_dry = 0.8d0
1646          aspect     = aspect_dry
1647          go to 1
1648       end if
1649
1650       if(fract_mass_water.lt.0.2d0) then
1651          aspect=aspect_dry-5.0d0*(aspect_dry-0.8d0)*fract_mass_water
1652          go to 1
1653       end if
1654
1655       if(fract_mass_water.ge.0.2d0.and.fract_mass_water.lt.0.8d0) &
1656       then
1657          aspect = 0.88d0-0.4d0*fract_mass_water
1658          goto 1
1659       end if
1660
1661       dd_melt = 1.d1*((6.0D0/pi)*hail_mass/den_water)**degree
1662
1663       aspect_melt = 0.9951d0+0.0251d0*dd_melt-0.03644*dd_melt**2+ &
1664                     0.005303*dd_melt**3-0.0002492*dd_melt**4
1665
1666       aspect      = 2.8d0-4.0d0*aspect_melt+5.0d0* &
1667                     (aspect_melt-0.56d0)*fract_mass_water
1668    1  continue
1669 !      dd=1.d1*(hail_mass*(1.0d0-fract_mass_water)/den0)**degree*factor
1670       fvw=den0*fract_mass_water/((1-fract_mass_water)*den_water+fract_mass_water*den0)
1671       dd=1.d1*((6.0D0/pi)*hail_mass/(fvw*den_water+(1.0d0-fvw)*den0))**degree
1672
1673       if(aspect.eq.1.0d0) aspect=0.9999d0
1674
1675       rp = dd*abs(sqrt(dc_wet))/(1.d1*wl)
1676
1677       if(rp.gt.0.1d0) then
1678          if(twolayer_hail == 1) then
1679              aspect2 = aspect
1680              dcore = (1.d0-fvw)**(1.d0/3.d0)*dd
1681              rpquad = rpquada(5)
1682              if ((rp.lt.rpquad) .OR. (usequad .EQV. .FALSE.)) then
1683             !call t_matrix2_dp(wl,dd,dcore,aspect,aspect2,dc_water,dc_hail,f_a,f_b,f_a0,f_b0)
1684                                     ! [KS] >> This code is not linked
1685              else
1686             !call t_matrix2_qp(wl,dd,dcore,aspect,aspect2,dc_water,dc_hail,f_a,f_b,f_a0,f_b0)
1687                                     ! [KS] >> This code is not linked
1688              endif
1689          else
1690              angle = 1.8d2
1691             ! call t_matrix(dd,wl,dc_wet,aspect,angle,f_a,f_b,6,'Hail')
1692              angle = 0.0d0
1693             ! call t_matrix(dd,wl,dc_wet,aspect,angle,f_a0,f_b0,6,'Hail')
1694              f_b0 = -DCONJG(f_b0)
1695              f_a0 = DCONJG(f_a0)
1696          endif
1697       else
1698 !     **************   external spheroid    ***************
1699
1700           dd2=1.d1*((6.0D0/pi)*hail_mass*(fract_mass_water/den_water+(1-fract_mass_water)/den0))**degree
1701
1702
1703           ff       = sqrt((1.0d0/aspect)**2-1.0d0)
1704           ff2      = ff**2
1705
1706           shape2_a = ((1.0+ff2)/ff2)*(1.0d0-atan(ff)/ff)
1707           shape2_b = 0.5d0*(1.0d0-shape2_a)
1708
1709 !     ************   inner spheroid   ********************
1710
1711           dd1=1.d1*((6.0D0/pi)*hail_mass*(1.0d0-fract_mass_water)/den0)**degree
1712
1713           shape1_a = shape2_a
1714           shape1_b = shape2_b
1715           psi      = (dd1/dd2)**3
1716
1717 !     ******   Scattering amplitudes    *********
1718
1719           tmp   =  pi**2*dd2**3/(6.0d2*wl**2)
1720
1721           num   = (dc_water-1.0d0)*(dc_water+(dc_hail-dc_water)* &
1722                   (shape1_a-psi*shape2_a))+ &
1723                   psi*dc_water*(dc_hail-dc_water)
1724
1725           denum = (dc_water+(dc_hail-dc_water)*(shape1_a-psi*shape2_a))* &
1726                   (1.0d0+(dc_water-1.0d0)*shape2_a)+ &
1727                   psi*shape2_a*dc_water*(dc_hail-dc_water)
1728
1729           f_a0   = tmp*num/denum
1730           f_a  = DCONJG(f_a0)
1731
1732           num   = (dc_water-1.0d0)*(dc_water+(dc_hail-dc_water)* &
1733                   (shape1_b-psi*shape2_b))+ &
1734                    psi*dc_water*(dc_hail-dc_water)
1735
1736           denum=(dc_water+(dc_hail-dc_water)*(shape1_b-psi*shape2_b))* &
1737                 (1.0d0+(dc_water-1.0d0)*shape2_b)+ &
1738                  psi*shape2_b*dc_water*(dc_hail-dc_water)
1739
1740           f_b0  = tmp*num/denum
1741           f_b = DCONJG(f_b0)
1742       end if
1743
1744       return
1745       end subroutine calc_scattering_hail
1746
1747 ! subroutine calc_scattering_hail
1748 !**** ********************************************************************
1749
1750       subroutine calc_scattering_fd &
1751                         (wl, fd_mass, den_bulk,fract_mass_water,dc_water,dc_ice,dc_wet, &
1752                          f_a,f_b,f_a0,f_b0)
1753
1754           USE scatt_tables,ONLY:twolayer_fd,rpquada,usequad
1755       !USE t_matrix2_quad_mod              ! ### [KS] : this is not linked since we use look-up-tables
1756       !USE t_matrix2_double_mod                ! ### [KS] : this is not linked since we use look-up-tables
1757
1758       !implicit double precision (a-h,o-z)
1759           implicit none
1760
1761           intrinsic DCONJG
1762
1763 ! ### Interface
1764           double precision :: wl, fd_mass, den_bulk, fract_mass_water
1765       complex(8) :: dc_water,dc_ice,dc_wet,num,denum,f_a,f_b,f_a0,f_b0
1766 ! ### INterfcae
1767
1768       double precision, parameter :: pi = 3.14159265D0, den_water = 1.0d0, den_ice = 0.91d0
1769           double precision :: degree, dd_dry, aspect_dry, aspect, aspect2, dd_melt, fvw, aspect_melt, dd, rp, &
1770                                                   angle, dd2, asp_w, asp_i, ff, ff2, shape2_a, shape2_b, dd1, shape1_a, shape1_b, &
1771                                                   psi, tmp, dcore, rpquad
1772
1773
1774       degree=1.0d0/3.0d0
1775
1776       dd_dry = 1.d1*((6.0D0/pi)*fd_mass/den_ice)**degree
1777
1778       if(dd_dry.lt.1.0d1) then
1779          aspect_dry = 1.0d0-2.0d-2*dd_dry
1780       else
1781          aspect_dry = 0.8d0
1782          aspect     = aspect_dry
1783          go to 1
1784       end if
1785
1786       if(fract_mass_water.lt.0.2d0) then
1787          aspect=aspect_dry-5.0d0*(aspect_dry-0.8d0)*fract_mass_water
1788          go to 1
1789       end if
1790
1791       if(fract_mass_water.ge.0.2d0.and.fract_mass_water.lt.0.8d0) &
1792       then
1793          aspect = 0.88d0-0.4d0*fract_mass_water
1794          goto 1
1795       end if
1796
1797       dd_melt = 1.d1*((6.0D0/pi)*fd_mass/den_water)**degree
1798
1799       aspect_melt = 0.9951d0+0.0251d0*dd_melt-0.03644*dd_melt**2+ &
1800                     0.005303*dd_melt**3-0.0002492*dd_melt**4
1801
1802       aspect      = 2.8d0-4.0d0*aspect_melt+5.0d0* &
1803                     (aspect_melt-0.56d0)*fract_mass_water
1804    1  continue
1805       fvw=den_ice*fract_mass_water/((1-fract_mass_water)*den_water+fract_mass_water*den_ice)
1806       dd=1.d1*((6.0D0/pi)*fd_mass/(fvw*den_water+(1.0d0-fvw)*den_ice))**degree
1807
1808       if(aspect.eq.1.0d0) aspect=0.9999d0
1809
1810       rp = dd*abs(sqrt(dc_wet))/(1.d1*wl)
1811
1812       if(rp.gt.0.1d0) then
1813          if(twolayer_fd == 1 ) then
1814              aspect2 = aspect
1815              dcore = fvw**(1.d0/3.d0)*dd
1816              rpquad = rpquada(2)
1817              if ((rp.lt.rpquad) .OR. (usequad .EQV. .FALSE.)) then
1818              !call t_matrix2_dp(wl,dd,dcore,aspect,aspect2,dc_ice,dc_water,f_a,f_b,f_a0,f_b0)
1819              ! [KS] >> This part is not reached
1820                          else
1821             !call t_matrix2_qp(wl,dd,dcore,aspect,aspect2,dc_ice,dc_water,f_a,f_b,f_a0,f_b0)
1822             ! >>[KS] : This part is not reached
1823                          endif
1824          else
1825              angle = 1.8d2
1826           !   call t_matrix(dd,wl,dc_wet,aspect,angle,f_a,f_b,4,'FD')
1827              angle = 0.0d0
1828           !   call t_matrix(dd,wl,dc_wet,aspect,angle,f_a0,f_b0,4,'FD')
1829              f_b0 = -DCONJG(f_b0)
1830              f_a0 = DCONJG(f_a0)
1831          endif
1832       else
1833
1834 !     **************   external spheroid    ***************
1835           dd2=1.d1*((6.0D0/pi)*fd_mass*(fract_mass_water/den_water+ &
1836              (1.0d0-fract_mass_water)/den_ice))**degree
1837
1838           if(dd2.lt.1.0d1) then
1839             asp_w = 0.9951d0+0.0251d0*dd2-0.03644*dd2**2+ &
1840                     0.005303*dd2**3-0.0002492*dd2**4
1841             asp_i = 1.0d0-0.02d0*dd2
1842           else
1843             asp_w = 0.4131d0
1844             asp_i = 0.8d0
1845           endif
1846
1847           aspect = fract_mass_water*asp_w+ &
1848                    (1.d0-fract_mass_water)*asp_i
1849
1850           if(aspect.eq.1.0d0) aspect=0.9999d0
1851
1852
1853           ff       = sqrt((1.0d0/aspect)**2-1.0d0)
1854           ff2      = ff**2
1855
1856           shape2_a = ((1.0+ff2)/ff2)*(1.0d0-atan(ff)/ff)
1857           shape2_b = 0.5d0*(1.0d0-shape2_a)
1858
1859
1860     !     ************   inner spheroid   ********************
1861
1862
1863           dd1=1.d1*((6.0D0/pi)*fd_mass*fract_mass_water/den_water)**degree
1864
1865           shape1_a = shape2_a
1866           shape1_b = shape2_b
1867           psi      = (dd1/dd2)**3
1868
1869
1870     !     ******   Scattering amplitudes    *********
1871
1872           tmp   =  pi**2*dd2**3/(6.0d2*wl**2)
1873
1874
1875           num   = (dc_ice-1.0d0)*(dc_ice+(dc_water-dc_ice)* &
1876                   (shape1_a-psi*shape2_a))+ &
1877                    psi*dc_ice*(dc_water-dc_ice)
1878
1879           denum = (dc_ice+(dc_water-dc_ice)*(shape1_a-psi*shape2_a))* &
1880                   (1.0d0+(dc_ice-1.0d0)*shape2_a)+ &
1881                    psi*shape2_a*dc_ice*(dc_water-dc_ice)
1882
1883
1884
1885           f_a0   = tmp*num/denum
1886           f_a  = DCONJG(f_a0)
1887
1888
1889           num   = (dc_ice-1.0d0)*(dc_ice+(dc_water-dc_ice)* &
1890                   (shape1_b-psi*shape2_b))+ &
1891                    psi*dc_ice*(dc_water-dc_ice)
1892
1893           denum=(dc_ice+(dc_water-dc_ice)*(shape1_b-psi*shape2_b))* &
1894                 (1.0d0+(dc_ice-1.0d0)*shape2_b)+ &
1895                  psi*shape2_b*dc_ice*(dc_water-dc_ice)
1896
1897           f_b0  = tmp*num/denum
1898           f_b = DCONJG(f_b0)
1899       endif
1900
1901       return
1902       end  subroutine calc_scattering_fd
1903 !**************************************************************************
1904       subroutine calc_scattering_grau1 &
1905                  (wl,grau_mass,den_bulk,fract_mass_water,dc,f_a,f_b,f_a0,f_b0)
1906
1907       !implicit double precision (a-h,o-z)
1908            implicit none
1909
1910            intrinsic DCONJG
1911 ! ### Interface
1912       double precision :: wl, grau_mass,den_bulk, fract_mass_water
1913           complex(8) :: dc, f_a, f_b, f_a0, f_b0
1914 ! ### Interface
1915
1916       double precision, parameter :: pi = 3.14159265D0, den0 = 0.4d0
1917       double precision,parameter :: den_water = 1.0d0
1918           double precision :: degree, dd, dd_dry, aspect_dry, aspect_melt, rp, angle, ff, ff2, shape_a, &
1919                                                   shape_b, tmp, aspect, dd_melt
1920
1921       degree=1.0d0/3.0d0
1922
1923      !dd=1.d1*(grau_mass*(1.0d0-fract_mass_water)/den0)**degree*factor
1924       dd=1.d1*((6.0D0/pi)*grau_mass*((1.0d0-fract_mass_water)/den0+fract_mass_water/den_water))**degree
1925
1926       dd_dry = 1.d1*((6.0D0/pi)*grau_mass/den0)**degree
1927
1928       if(dd_dry.lt.1.0d1) then
1929          aspect_dry = 1.0d0-2.0d-2*dd_dry
1930       else
1931          aspect_dry = 0.8d0
1932       endif
1933
1934       if(fract_mass_water.lt.0.2d0) then
1935          aspect=aspect_dry-5.0d0*(aspect_dry-0.8d0)*fract_mass_water
1936          goto 1
1937       endif
1938
1939       if(fract_mass_water.ge.0.2d0.and.fract_mass_water.lt.0.8d0) &
1940       then
1941          aspect = 0.88d0-0.4d0*fract_mass_water
1942          goto 1
1943       endif
1944
1945       dd_melt     = 1.d1*((6.0D0/pi)*grau_mass/den_water)**degree
1946
1947       if(dd_melt.lt.1.0d1) then
1948
1949          aspect_melt=0.9951d0+0.0251d0*dd_melt-0.03644*dd_melt**2+ &
1950                      0.005303*dd_melt**3-0.0002492*dd_melt**4
1951       else
1952          aspect_melt  = 0.4131d0
1953       endif
1954
1955       aspect = 2.8d0-4.0d0*aspect_melt+5.0d0* &
1956               (aspect_melt-0.56d0)*fract_mass_water
1957    1  continue
1958
1959       if(aspect.eq.1.0d0) aspect=0.9999d0
1960
1961       rp = dd*abs(sqrt(dc))/(1.d1*wl)
1962
1963       if(rp.gt.0.1d0) then
1964
1965          angle = 1.8d2
1966
1967         ! call t_matrix(dd,wl,dc,aspect,angle,f_a,f_b,7,'grau1') ! [KS] >> This is not linked as we use lookup tables
1968
1969          angle = 0.0d0
1970
1971          ! call t_matrix(dd,wl,dc,aspect,angle,f_a0,f_b0,7,'grau1') ! [KS] >> This is not linked as we use lookup tables
1972          f_b0 = -DCONJG(f_b0)
1973          f_a0 = DCONJG(f_a0)
1974
1975       else
1976
1977          ff      = sqrt((1.0d0/aspect)**2-1.0d0)
1978          ff2     = ff**2
1979          shape_a = ((1.0+ff2)/ff2)*(1.0d0-atan(ff)/ff)
1980          shape_b = 0.5d0*(1.0d0-shape_a)
1981          tmp     = pi**2*dd**3/(6.0d2*wl**2)
1982          f_a0    = tmp/(shape_a+1.0d0/(dc-1.0d0))
1983          f_b0    = tmp/(shape_b+1.0d0/(dc-1.0d0))
1984          f_a     = DCONJG(f_a0)
1985          f_b     = DCONJG(f_b0)
1986
1987       end if
1988
1989       return
1990       end subroutine calc_scattering_grau1
1991
1992 ! subroutine calc_scattering_grau1
1993 !**** ********************************************************************
1994 ! JCS -- Modified list of arguments and equation for dd.  Now, dc_wet_inner
1995 ! represents the dc of the inner spongy 'layer', whereas dc_wet represents the
1996 ! dc for the entire particle. We need both because the t_matrix scattering
1997 ! subroutine still assume a homogeneous mixture; in contast, the two-layer
1998 ! Rayleigh calculations work on the spongy inner layer, for which dc_wet_inner
1999 ! is used.
2000       subroutine calc_scattering_grau2 &
2001                                 (wl,grau_mass,den_inner, &
2002                                              fract_mass_water,fract_mass_crit,dc_water,dc_ice,dc_wet_inner,dc_wet, &
2003                                  f_a,f_b,f_a0,f_b0)
2004
2005       USE scatt_tables,ONLY:twolayer_graupel,rpquada,usequad
2006       !USE t_matrix2_quad_mod                ! ### [KS] : this is not linked since we use look-up-tables
2007       !USE t_matrix2_double_mod              ! ### [KS] : this is not linked since we use look-up-tables
2008
2009           !implicit double precision (a-h,o-z)
2010           implicit none
2011
2012           intrinsic DCONJG
2013 ! ### INterface
2014           double precision :: wl,grau_mass,den_inner, fract_mass_water,fract_mass_crit
2015       complex(8) :: dc_water,dc_ice,dc_wet_inner,dc_wet,num,denum,f_a,f_b,f_a0,f_b0
2016 ! ### Interface
2017       double precision,parameter :: pi = 3.14159265D0, den_water = 1.0d0, den0 = 0.4d0
2018           double precision :: degree, dd, dd_dry, aspect_dry, aspect, aspect2, dd_melt, rp, angle, &
2019                                              fract_mass_excess, shape2_a, shape2_b, dd1, psi, tmp, aspect_melt, &
2020                                                  dd2, ff, ff2, shape1_a, shape1_b, rpquad, fvw
2021
2022       degree=1.0d0/3.0d0
2023
2024 !      dd=1.d1*(grau_mass*(1.0d0-fract_mass_water)/den0)**degree*factor
2025       dd=1.d1*((6.0D0/pi)*grau_mass*((1.0d0-fract_mass_water)/den0+fract_mass_water/den_water))**degree
2026
2027       dd_dry = 1.d1*((6.0D0/pi)*grau_mass/den0)**degree
2028
2029       if(dd_dry.lt.1.0d1) then
2030          aspect_dry = 1.0d0-2.0d-2*dd_dry
2031       else
2032          aspect_dry = 0.8d0
2033       endif
2034
2035       if(fract_mass_water.lt.0.2d0) then
2036          aspect=aspect_dry-5.0d0*(aspect_dry-0.8d0)*fract_mass_water
2037          goto 1
2038       endif
2039
2040       if(fract_mass_water.ge.0.2d0.and.fract_mass_water.lt.0.8d0) &
2041       then
2042          aspect = 0.88d0-0.4d0*fract_mass_water
2043          goto 1
2044       endif
2045
2046       dd_melt     = 1.d1*((6.0D0/pi)*grau_mass/den_water)**degree
2047
2048       if(dd_melt.lt.1.0d1) then
2049
2050          aspect_melt=0.9951d0+0.0251d0*dd_melt-0.03644*dd_melt**2+ &
2051                      0.005303*dd_melt**3-0.0002492*dd_melt**4
2052       else
2053          aspect_melt  = 0.4131d0
2054       endif
2055
2056       aspect = 2.8d0-4.0d0*aspect_melt+5.0d0* &
2057               (aspect_melt-0.56d0)*fract_mass_water
2058    1  continue
2059 ! JCS - fract_mass_excess is the excess water that will be kept entirely water.
2060 ! The rest of the water was used to soak the interior of the particle.
2061       fract_mass_excess=fract_mass_water-fract_mass_crit
2062 ! dd2 is outer diameter
2063       dd2=1.d1*((6.0D0/pi)*grau_mass*(fract_mass_excess/den_water+ &
2064          (1.0d0-fract_mass_excess)/den_inner))**degree
2065 ! dd1 is core or inner diameter
2066       dd1=1.d1*((6.0D0/pi)*grau_mass*(1.0d0-fract_mass_excess)/den_inner)**degree
2067
2068       if(aspect.eq.1.0d0) aspect=0.9999d0
2069
2070       rp = dd*abs(sqrt(dc_wet))/(1.d1*wl)
2071
2072       if(rp.gt.0.1d0) then
2073          if(twolayer_graupel == 1) then
2074              aspect2 = aspect
2075              fvw = den_inner*fract_mass_excess/((1.d0-fract_mass_excess)*den_water+fract_mass_excess*den_inner)
2076              rpquad = rpquada(4)
2077              if ((rp.lt.rpquad) .OR. (usequad .EQV. .FALSE.)) then
2078                  !call t_matrix2_dp(wl,dd2,dd1,aspect,aspect2,dc_water,dc_wet_inner,f_a,f_b,f_a0,f_b0)
2079                    ! [KS] >> This part is not reached
2080                          else
2081                  !call t_matrix2_qp(wl,dd2,dd1,aspect,aspect2,dc_water,dc_wet_inner,f_a,f_b,f_a0,f_b0)
2082                   ! [KS] >> This part is not reached
2083                          endif
2084          else
2085              angle = 1.8d2
2086             ! call t_matrix(dd,wl,dc_wet,aspect,angle,f_a,f_b,7,'grau2') ! [KS] >> This is not linked as we use lookup tables
2087              angle = 0.0d0
2088             ! call t_matrix(dd,wl,dc_wet,aspect,angle,f_a0,f_b0,7,'grau2') ! [KS] >> This is not linked as we use lookup tables
2089              f_b0 = -DCONJG(f_b0)
2090              f_a0 = DCONJG(f_a0)
2091          endif
2092       else
2093
2094         !    **************   external spheroid    ***************
2095
2096           ff  = sqrt((1.0d0/aspect)**2-1.0d0)
2097           ff2 = ff**2
2098
2099           shape2_a = ((1.0+ff2)/ff2)*(1.0d0-atan(ff)/ff)
2100
2101           shape2_b = 0.5d0*(1.0d0-shape2_a)
2102
2103     !     ************   Inner spheroid    ***************
2104
2105           shape1_a =  shape2_a
2106           shape1_b =  shape2_b
2107
2108           psi      = (dd1/dd2)**3
2109
2110     !      ********  Scettering amplitudes    ****
2111
2112           tmp   = pi**2*dd2**3/(6.0d2*wl**2)
2113
2114           num   = (dc_water-1.0d0)*(dc_water+(dc_wet_inner-dc_water)* &
2115                   (shape1_a-psi*shape2_a))+ &
2116                    psi*dc_water*(dc_wet_inner-dc_water)
2117
2118           denum = (dc_water+(dc_wet_inner-dc_water)* &
2119                   (shape1_a-psi*shape2_a))* &
2120                   (1.0d0+(dc_water-1.0d0)*shape2_a)+ &
2121                    psi*shape2_a*dc_water*(dc_wet_inner-dc_water)
2122
2123           f_a0   = tmp*num/denum
2124           f_a  = DCONJG(f_a0)
2125
2126           num   = (dc_water-1.0d0)*(dc_water+(dc_wet_inner-dc_water)* &
2127                   (shape1_b-psi*shape2_b))+ &
2128                    psi*dc_water*(dc_wet_inner-dc_water)
2129
2130           denum = (dc_water+(dc_wet_inner-dc_water)* &
2131                   (shape1_b-psi*shape2_b))* &
2132                   (1.0d0+(dc_water-1.0d0)*shape2_b)+ &
2133                    psi*shape2_b*dc_water*(dc_wet_inner-dc_water)
2134
2135           f_b0   = tmp*num/denum
2136           f_b  = DCONJG(f_b0)
2137       endif
2138       return
2139       end subroutine calc_scattering_grau2
2140
2141 ! subroutine calc_scattering_grau2
2142 !**** ********************************************************************
2143       subroutine calc_rayleigh_plate &
2144                                                 (wl,plate_mass,den_bulk, &
2145                                              fract_mass_water,dc_water,dc_plate,f_a,f_b,f_a0,f_b0)
2146
2147       ! implicit double precision (a-h,o-z)
2148                 implicit none
2149
2150                 intrinsic DCONJG
2151 ! ### Interface
2152       double precision :: wl,plate_mass,den_bulk, fract_mass_water
2153           complex(8) :: dc_water,dc_plate,num,denum,f_a,f_b,f_a0,f_b0
2154 ! ### INterface
2155
2156       double precision, parameter :: pi = 3.14159265D0, den_water = 1.0d0
2157       double precision, parameter :: alfa  = 0.047d0, beta = 0.474d0
2158           double precision :: degree, dd2, dd_dry, tmp, aa, bb, aspect_dry, dd_melt, aspect, aspect_melt, &
2159                                                   ff, ff2, shape2_a, shape2_b, psi, dd1, shape1_a, shape1_b
2160
2161       degree=1.0d0/3.0d0
2162
2163 !     **************   external spheroid    ***************
2164
2165       dd2=1.d1*((6.0D0/pi)*plate_mass*(fract_mass_water/den_water+ &
2166          (1.0d0-fract_mass_water)/den_bulk))**degree
2167
2168       dd_dry     = 1.d1*((6.0D0/pi)*plate_mass/den_bulk)**degree
2169       tmp        = dd_dry**3/alfa
2170       aa         = alfa*tmp**(beta/(2.0d0+beta))
2171       bb         = tmp**(1.0d0/(2.0d0+beta))
2172       aspect_dry = aa/bb
2173
2174       if(aspect_dry.gt.1.0d0) aspect_dry = 1.0d0
2175
2176       dd_melt = 1.d1*((6.0D0/pi)*plate_mass/den_water)**degree
2177
2178       if(dd_melt.lt.1.0d1) then
2179          aspect_melt=0.9951d0+0.0251d0*dd_melt-0.03644*dd_melt**2+ &
2180                      0.005303*dd_melt**3-0.0002492*dd_melt**4
2181       else
2182          aspect_melt  = 0.4131d0
2183       endif
2184
2185       aspect=aspect_dry+fract_mass_water*(aspect_melt-aspect_dry)
2186
2187       if(aspect.eq.1.0d0) aspect=0.9999d0
2188
2189       ff  = sqrt((1.0d0/aspect)**2-1.0d0)
2190       ff2 = ff**2
2191
2192       shape2_a = ((1.0+ff2)/ff2)*(1.0d0-atan(ff)/ff)
2193
2194       shape2_b = 0.5d0*(1.0d0-shape2_a)
2195
2196 !     ************   inner spheroid   ********************
2197
2198       dd1= &
2199       1.d1*((6.0D0/pi)*plate_mass*(1.0d0-fract_mass_water)/den_bulk)**degree
2200
2201       shape1_a =  shape2_a
2202
2203       shape1_b =  shape2_b
2204
2205       psi      = (dd1/dd2)**3
2206
2207 !     ******   Scattering amplitudes    *********
2208
2209       tmp = pi**2*dd2**3/(6.0d2*wl**2)
2210
2211       num = (dc_water-1.0d0)*(dc_water+(dc_plate-dc_water)* &
2212             (shape1_a-psi*shape2_a))+ &
2213              psi*dc_water*(dc_plate-dc_water)
2214
2215       denum=(dc_water+(dc_plate-dc_water)*(shape1_a-psi*shape2_a))* &
2216             (1.0d0+(dc_water-1.0d0)*shape2_a)+ &
2217              psi*shape2_a*dc_water*(dc_plate-dc_water)
2218
2219       f_a0   = tmp*num/denum
2220       f_a  = DCONJG(f_a0)
2221
2222       num = (dc_water-1.0d0)*(dc_water+(dc_plate-dc_water)* &
2223             (shape1_b-psi*shape2_b))+ &
2224              psi*dc_water*(dc_plate-dc_water)
2225
2226       denum=(dc_water+(dc_plate-dc_water)*(shape1_b-psi*shape2_b))* &
2227             (1.0d0+(dc_water-1.0d0)*shape2_b)+ &
2228              psi*shape2_b*dc_water*(dc_plate-dc_water)
2229
2230       f_b0   = tmp*num/denum
2231       f_b  = DCONJG(f_b0)
2232
2233       return
2234       end subroutine calc_rayleigh_plate
2235
2236 ! subroutine calc_rayleigh_plate &
2237
2238 !**** *************************************************************** &
2239 !**** ***************************************************************
2240
2241       subroutine calc_rayleigh_dendr &
2242                                                 (wl,dendr_mass,den_bulk, &
2243                                                  fract_mass_water,dd_dry,dc,f_a,f_b,f_a0,f_b0,ijk,kx,kz,kb)
2244
2245       !implicit double precision (a-h,o-z)
2246            implicit none
2247
2248            intrinsic DCONJG
2249 ! ### Interface
2250           integer :: ijk, kx, kz, kb
2251           double precision :: wl,dendr_mass,den_bulk, fract_mass_water,dd_dry
2252       complex(8) :: dc, f_a, f_b, f_a0, f_b0
2253 ! ### Interface
2254       double precision,parameter :: pi = 3.14159265D0, den_water = 1.0d0
2255       double precision,parameter :: alfa = 0.038d0, beta = 0.377d0
2256           double precision :: degree, dd, tmp, aa, bb, aspect_dry, dD_melt, aspect_melt, aspect, &
2257                                                   rp, angle, ff, ff2, shpae_a, shape_b, c, shape_a
2258
2259       degree=1.0d0/3.0d0
2260
2261       dd = 1.d1*((6.0D0/pi)*dendr_mass*(fract_mass_water/den_water+ &
2262           (1.0d0-fract_mass_water)/den_bulk))**degree
2263
2264       tmp         = dd_dry**3/alfa
2265       aa          = alfa*tmp**(beta/(2.0d0+beta))
2266       bb          = tmp**(1.0d0/(2.0d0+beta))
2267       aspect_dry  = aa/bb
2268
2269       if(aspect_dry.gt.1.0d0) aspect_dry = 1.0d0
2270
2271       dd_melt = 1.d1*((6.0D0/pi)*dendr_mass/den_water)**degree
2272
2273       if(dd_melt.lt.1.0d1) then
2274          aspect_melt=0.9951d0+0.0251d0*dd_melt-0.03644*dd_melt**2+ &
2275                      0.005303*dd_melt**3-0.0002492*dd_melt**4
2276       else
2277          aspect_melt  = 0.4131d0
2278       end if
2279
2280       aspect  = aspect_dry+fract_mass_water*(aspect_melt-aspect_dry)
2281
2282       if(aspect.eq.1.0d0) aspect=0.9999d0
2283
2284       rp = dd*abs(sqrt(dc))/(1.d1*wl)
2285
2286       if(rp.gt.0.4d0) then
2287
2288          angle = 1.8d2
2289
2290         ! call t_matrix(dd, wl,dc,aspect,angle,f_a,f_b,7,'dendr') ! [KS] >> This is not linked as we use lookup tables
2291
2292          angle = 0.0d0
2293
2294          ! call t_matrix(dd,wl,dc,aspect,angle,f_a0,f_b0,7,'dendr') ! [KS] >> This is not linked as we use lookup tables
2295          f_b0 = -DCONJG(f_b0)
2296          f_a0 = DCONJG(f_a0)
2297
2298 ! in case rp.gt.0.4d0
2299
2300       else
2301
2302 ! in case rp.le.0.4d0
2303
2304          ff      = sqrt((1.0d0/aspect)**2-1.0d0)
2305          ff2     = ff**2
2306
2307 ! new change 11.07.08                                         (start)
2308              c=atan(ff)/ff
2309 ! new change 11.07.08                                           (end)
2310
2311          shape_a = ((1.0d0+ff2)/ff2)*(1.0d0-atan(ff)/ff)
2312
2313          shape_b = 0.5d0*(1.0d0-shape_a)
2314
2315          tmp     = pi**2*dd**3/(6.0d2*wl**2)
2316
2317          f_a0     = tmp/(shape_a+1.0d0/(dc-1.0d0))
2318          f_a    = DCONJG(f_a0)
2319          f_b0     = tmp/(shape_b+1.0d0/(dc-1.0d0))
2320          f_b    = DCONJG(f_b0)
2321
2322       endif
2323
2324       return
2325       end subroutine calc_rayleigh_dendr
2326
2327 ! subroutine calc_rayleigh_dendr &
2328
2329 !**** ******************************************************************** &
2330 !**** ********************************************************************
2331
2332       subroutine calc_rayleigh_snow &
2333                         (wl,snow_mass,den_bulk, &
2334                          fract_mass_water,dd_dry,dc,f_a,f_b,f_a0,f_b0)
2335
2336       !implicit double precision (a-h,o-z)
2337           implicit none
2338
2339           intrinsic DCONJG
2340 ! ### INterface
2341       double precision :: wl,snow_mass,den_bulk, fract_mass_water,dd_dry
2342       complex(8) ::       dc, f_a, f_b, f_a0, f_b0
2343 ! ### INterface
2344       double precision, parameter :: pi = 3.14159265D0, den_water = 1.0d0
2345           double precision :: degree, dd, aspect_dry, dD_melt, aspect_melt, aspect, rp, angle, ff, ff2, &
2346                                                   shape_a, shape_b, tmp
2347
2348       degree=1.0d0/3.0d0
2349
2350       dd=1.d1*((6.0D0/pi)*snow_mass*(fract_mass_water/den_water+ &
2351         (1.0d0-fract_mass_water)/den_bulk))**degree
2352
2353       if(dd_dry.lt.1.0d1) then
2354          aspect_dry = 1.0d0-2.0d-2*dd_dry
2355       else
2356          aspect_dry = 0.8d0
2357       end if
2358
2359       dd_melt = 1.d1*((6.0D0/pi)*snow_mass/den_water)**degree
2360
2361       if(dd_melt.lt.1.0d1) then
2362          aspect_melt=0.9951d0+0.0251d0*dd_melt-0.03644*dd_melt**2+ &
2363                      0.005303*dd_melt**3-0.0002492*dd_melt**4
2364       else
2365          aspect_melt  = 0.4131d0
2366       end if
2367
2368       aspect = aspect_dry+fract_mass_water*(aspect_melt-aspect_dry)
2369
2370       if(aspect.eq.1.0d0) aspect=0.9999d0
2371
2372       rp = dd*abs(sqrt(dc))/(1.d1*wl)
2373
2374       if(rp.gt.0.1d0) then
2375
2376          angle = 1.8d2
2377
2378         ! call t_matrix(dd,wl,dc,aspect,angle,f_a,f_b,6,'snow') ! [KS] >> This is not linked as we use lookup tables
2379
2380          angle = 0.0d0
2381
2382          ! call t_matrix(dd,wl,dc,aspect,angle,f_a0,f_b0,6,'snow') ! [KS] >> This is not linked as we use lookup tables
2383          f_b0 = -DCONJG(f_b0)
2384          f_a0 = DCONJG(f_a0)
2385       else
2386
2387          ff      = sqrt((1.0d0/aspect)**2-1.0d0)
2388          ff2     = ff**2
2389
2390          shape_a = ((1.0+ff2)/ff2)*(1.0d0-atan(ff)/ff)
2391
2392          shape_b = 0.5d0*(1.0d0-shape_a)
2393
2394          tmp     = pi**2*dd**3/(6.0d2*wl**2)
2395
2396          f_a0     = tmp/(shape_a+1.0d0/(dc-1.0d0))
2397          f_b0     = tmp/(shape_b+1.0d0/(dc-1.0d0))
2398          f_a    = DCONJG(f_a0)
2399          f_b    = DCONJG(f_b0)
2400
2401       endif
2402
2403       return
2404       end subroutine calc_rayleigh_snow
2405
2406 ! subroutine calc_rayleigh_snow &
2407
2408 !**** ******************************************************************** &
2409 !**** ********************************************************************
2410
2411       subroutine calc_rayleigh_colum &
2412                         (wl,colum_mass,den_bulk, &
2413                         fract_mass_water,dd_dry,dc_water,dc_colum, &
2414                         f_a,f_b,f_a0,f_b0)
2415
2416       !implicit double precision (a-h,o-z)
2417            implicit none
2418
2419            intrinsic DCONJG
2420 ! ### Interface
2421           double precision :: wl,colum_mass,den_bulk, fract_mass_water,dd_dry
2422       complex(8) :: dc_water, dc_colum, num, denum, f_a, f_b, f_a0, f_b0
2423 ! ### Interface
2424       double precision,parameter :: pi  = 3.14159265D0, den_water = 1.0d0
2425       double precision, parameter :: alfa = 0.308d0, beta = 0.927d0
2426           double precision :: degree, dd2, tmp, aa, bb, aspect_dry, dd_melt, aspect_melt, ff, ff2, &
2427                                                   shape2_a, shape2_b, psi, aspect, dd1, shape1_a, shape1_b
2428
2429
2430       degree=1.0d0/3.0d0
2431 !     **************   external spheroid    ***************
2432
2433       dd2 = 1.d1*((6.0D0/pi)*colum_mass*(fract_mass_water/den_water+ &
2434            (1.0d0-fract_mass_water)/den_bulk))**degree
2435
2436       tmp        = dd_dry**3/alfa**2
2437
2438       aa         = tmp**(1.0d0/(2.0d0*beta+1.0d0))
2439
2440       bb         = alfa*tmp**(beta/(2.0d0*beta+1.0d0))
2441
2442       aspect_dry = aa/bb
2443
2444       dd_melt = 1.d1*((6.0D0/pi)*colum_mass/den_water)**degree
2445
2446       if(dd_melt.lt.1.0d1) then
2447          aspect_melt=0.9951d0+0.0251d0*dd_melt-0.03644*dd_melt**2+ &
2448                      0.005303*dd_melt**3-0.0002492*dd_melt**4
2449       else
2450          aspect_melt  = 0.4131d0
2451       endif
2452
2453       aspect = aspect_dry+fract_mass_water*(aspect_melt-aspect_dry)
2454
2455       if(aspect.eq.1.0d0) aspect=0.9999d0
2456
2457       if(aspect.lt.1.0d0) then
2458
2459          ff       = sqrt((1.0d0/aspect)**2-1.0d0)
2460          ff2      = ff**2
2461          shape2_a = ((1.0+ff2)/ff2)*(1.0d0-atan(ff)/ff)
2462          shape2_b = 0.5d0*(1.0d0-shape2_a)
2463
2464       else
2465
2466          ff       = sqrt(1.0d0-(1.0d0/aspect)**2)
2467          ff2      = ff**2
2468          shape2_a = ((1.0-ff2)/ff2)*(0.5d0* &
2469                     dlog((1.0d0+ff)/(1.0d0-ff))/ff-1.0d0)
2470          shape2_b = 0.5d0*(1.0d0-shape2_a)
2471
2472       endif
2473
2474 !     ************   inner spheroid   ********************
2475
2476       dd1=1.d1* &
2477       ((6.0D0/pi)*colum_mass*(1.0d0-fract_mass_water)/den_bulk)**degree
2478
2479       shape1_a = shape2_a
2480       shape1_b = shape2_b
2481       psi      = (dd1/dd2)**3
2482
2483 !     ******   Scattering amplitudes    *********
2484
2485       tmp   = pi**2*dd2**3/(6.0d2*wl**2)
2486
2487       num   = (dc_water-1.0d0)*(dc_water+(dc_colum-dc_water)* &
2488               (shape1_a-psi*shape2_a))+ &
2489               psi*dc_water*(dc_colum-dc_water)
2490
2491       denum = (dc_water+(dc_colum-dc_water)* &
2492               (shape1_a-psi*shape2_a))* &
2493               (1.0d0+(dc_water-1.0d0)*shape2_a)+ &
2494                psi*shape2_a*dc_water*(dc_colum-dc_water)
2495
2496       f_a0   = tmp*num/denum
2497       f_a  = DCONJG(f_a0)
2498
2499       num   = (dc_water-1.0d0)*(dc_water+(dc_colum-dc_water)* &
2500               (shape1_b-psi*shape2_b))+ &
2501               psi*dc_water*(dc_colum-dc_water)
2502
2503       denum = (dc_water+(dc_colum-dc_water)* &
2504               (shape1_b-psi*shape2_b))* &
2505               (1.0d0+(dc_water-1.0d0)*shape2_b)+ &
2506               psi*shape2_b*dc_water*(dc_colum-dc_water)
2507
2508       f_b0   = tmp*num/denum
2509       f_b  = DCONJG(f_b0)
2510
2511       return
2512       end subroutine calc_rayleigh_colum
2513
2514 ! subroutine calc_rayleigh_colum &
2515
2516 !**** ******************************************************************** &
2517 !**** ********************************************************************
2518
2519       subroutine calc_orient_colum(beta,a)
2520
2521       !implicit double precision (a-h,o-z)
2522            implicit none
2523
2524         double precision :: a(7), beta
2525
2526       a(1)=  0.5d0*sin(beta)**2
2527       a(2)=  0.5d0
2528       a(3)=  0.375d0*sin(beta)**2
2529       a(4)=  0.375d0
2530       a(5)=  0.125d0*sin(beta)**2
2531       a(6)=  0.0d0
2532       a(7)= -0.5d0*cos(beta)**2
2533
2534       return
2535       end subroutine calc_orient_colum
2536
2537 ! subroutine calc_orient_colum &
2538
2539 !**** ******************************************************************** &
2540 !**** ********************************************************************
2541 ! JCS - fixed a(2) consistent with Rhzykov et al. (2011)
2542       subroutine calc_orient_water(a)
2543
2544       !implicit double precision (a-h,o-z)
2545          implicit none
2546 ! ### Interface
2547          double precision :: a(7)
2548 ! ### Interface
2549
2550       !dimension a(7)
2551       double precision,parameter :: sigma_r = 0.17453292d0
2552           double precision :: sigma, r
2553
2554
2555       sigma   = sigma_r
2556       r       = dexp(-2.0d0*sigma**2)
2557
2558       a(1)= 0.25d0*(1.0+r)**2
2559       a(2)= 0.25d0*(1.0-r**2)
2560       a(3)= (0.375d0+0.5d0*r+0.125d0*r**4)**2
2561       a(4)= (0.375d0-0.5d0*r+0.125d0*r**4)* &
2562             (0.375d0+0.5d0*r+0.125d0*r**4)
2563       a(5)= 0.125d0*(0.375d0+0.5d0*r+0.125d0*r**4)*(1.0d0-r**4)
2564       a(6)= 0.0d0
2565       a(7)= 0.5d0*r*(1.0d0+r)
2566
2567       return
2568       end subroutine calc_orient_water
2569
2570 ! subroutine calc_orient_water &
2571
2572 !**** *************************************************************** &
2573 !**** ***************************************************************
2574 ! new change 4.08.11                                          (start)
2575 ! JCS - fixed a(2,kb) consistent with Rhzykov et al. (2011)
2576       subroutine calc_orient(fract_mass_water,a,kb,number_bin) ! ### (KS)
2577
2578 !use microprm
2579
2580       !implicit double precision (a-h,o-z)
2581           implicit none
2582 ! ### Interface
2583           integer :: number_bin, kb
2584       double precision :: a(7,number_bin), fract_mass_water
2585
2586 ! ### Interface
2587
2588       double precision,parameter :: sigma_r = 0.17453292d0, sigma_s  = 0.69813176d0
2589           double precision :: sigma, r
2590
2591
2592       sigma   = sigma_s+fract_mass_water*(sigma_r-sigma_s)
2593       r       = dexp(-2.0d0*sigma**2)
2594
2595       a(1,kb)= 0.25d0*(1.0+r)**2
2596       a(2,kb)= 0.25d0*(1.0-r**2)
2597       a(3,kb)=(0.375d0+0.5d0*r+0.125d0*r**4)**2
2598       a(4,kb)=(0.375d0-0.5d0*r+0.125d0*r**4)* &
2599               (0.375d0+0.5d0*r+0.125d0*r**4)
2600       a(5,kb)= 0.125d0*(0.375d0+0.5d0*r+0.125d0*r**4)*(1.0d0-r**4)
2601       a(6,kb)= 0.0d0
2602       a(7,kb)= 0.5d0*r*(1.0d0+r)
2603
2604       return
2605       end subroutine calc_orient
2606
2607 ! subroutine calc_orient
2608
2609 ! new change 4.08.11                                            (end) &
2610 !**** *************************************************************** &
2611 !**** ***************************************************************
2612 ! Andrei's new change of 4.08.11                              (start)
2613
2614       subroutine integr &
2615                                 (a,bin_conc,f_a,f_b,f_a0,f_b0,zh,zv,ldr,kdp,rhv,cdr,ah,adp,ijk, &
2616                                  kx,kz,ihydromet,number_bin)
2617
2618 !use microprm
2619
2620       !implicit double precision (a-h,o-z)
2621           implicit none
2622
2623           intrinsic dimag, dconjg
2624 ! ### Interface
2625                 ! parameter(number_bin = NKR_43Bins)
2626                 integer :: number_bin, ijk, kx, kz, ihydromet, kb
2627                 double precision :: ldr, kdp
2628                 complex(8) :: rhv
2629                 complex(8) :: f_a(number_bin),  f_b(number_bin), &
2630                                            f_a0(number_bin), f_b0(number_bin)
2631                 double precision :: a(7,number_bin)
2632                 double precision :: bin_conc(number_bin), zh, zv, cdr, ah, adp
2633 ! ### Interface
2634
2635                 double precision :: cdrn, cdrd, b
2636                 double precision :: aj1n, aj1d, aj1, aj2n, aj2d, aj2
2637
2638           zh  = 0.0d0
2639           zv  = 0.0d0
2640           ldr = 0.0d0
2641           kdp = 0.0d0
2642           rhv = (0.0d0,0.0d0)
2643           ah = 0.0d0
2644           adp = 0.0d0
2645           cdr = 0.0d0
2646           cdrn = 0.0d0
2647           cdrd = 0.0d0
2648           aj1n = 0.0d0
2649           aj1d = 0.0d0
2650           aj1 = 0.0d0
2651           aj2n = 0.d0
2652           aj2d = 0.d0
2653           aj2 = 0.d0
2654
2655           do kb=1,number_bin
2656
2657                         zh = zh + bin_conc(kb)*(abs(f_b(kb))**2- &
2658                         2.0d0*a(2,kb)*dble(dconjg(f_b(kb))*(f_b(kb)-f_a(kb)))+ &
2659                         a(4,kb)*abs(f_b(kb)-f_a(kb))**2)
2660
2661                 zv = zv + bin_conc(kb)*(abs(f_b(kb))**2- &
2662                 2.0d0*a(1,kb)*dble(dconjg(f_b(kb))*(f_b(kb)-f_a(kb)))+ &
2663                 a(3,kb)*abs(f_b(kb)-f_a(kb))**2)
2664
2665                 ldr = ldr + bin_conc(kb)*a(5,kb)*abs(f_b(kb)-f_a(kb))**2
2666
2667                         b = dble(f_b0(kb)-f_a0(kb))
2668
2669                 kdp = kdp + bin_conc(kb)*a(7,kb)*b
2670
2671                 rhv = rhv + bin_conc(kb)*(abs(f_b(kb))**2+ &
2672                         a(5,kb)*abs(f_b(kb)-f_a(kb))**2- &
2673                         a(1,kb)*dconjg(f_b(kb))*(f_b(kb)-f_a(kb))- &
2674                         a(2,kb)*f_b(kb)*dconjg(f_b(kb)-f_a(kb)))
2675
2676 ! JCS - CDR is from Eqn (17) in Ryzhkov (2001)
2677 ! THIS DOES NOT ACCOUNT FOR CANTING ANGLE VARIABILITY (ASSUMES VERY SMALL
2678 ! DISTRIBUTION OF CANTING ANGLES)
2679                 cdrn = cdrn + bin_conc(kb)*(abs(f_b(kb)-f_a(kb)))**2
2680                 cdrd = cdrd + bin_conc(kb)*(abs(f_b(kb)+f_a(kb)))**2
2681 ! JCS - Test -- CDR from Eqn (11) and (16) in Ryzhkov (2001) to account for
2682 ! canting angle variability
2683                         !aj1n = aj1n + bin_conc(kb)*sqrt(a(3,kb))*(abs(f_b(kb)-f_a(kb)))**2
2684                         !aj1d = aj1d + bin_conc(kb)*(abs(f_b(kb)))**2
2685                         !if (aj1d>0.d0) then
2686                         !    aj1=aj1+(aj1n/aj1d)
2687                         !endif
2688                         !aj2n = aj2n + bin_conc(kb)*(conjg(f_b(kb))*(f_b(kb)-f_a(kb)))
2689                         !aj2d = aj2d + bin_conc(kb)*(abs(f_a(kb)))**2
2690                         !if (aj2d>0.d0) then
2691                         !    aj2=aj2+(aj2n/aj2d)
2692                         !endif
2693 ! JCS - Equations for ah and adp are from Eqn 6 in Ryzhkov et al. (2013a)
2694          ah = ah + bin_conc(kb)*(dimag(f_b0(kb))-a(2,kb)*dimag(f_b0(kb)-f_a0(kb)))
2695
2696 ! JCS - Note that A5 in their paper is really a(7,kb) as defined in this program
2697          adp = adp + bin_conc(kb)*dimag(f_b0(kb)-f_a0(kb))*a(7,kb)
2698       enddo
2699
2700                 ! cdrn = (1.d0/4.d0)*aj1
2701                 ! cdrd = 1-dble(aj2)+1.d0/4.d0*(aj1)
2702         if (abs(cdrd)>0.d0) cdr = cdrn/cdrd
2703
2704           return
2705       end subroutine integr
2706
2707 ! subroutine integr
2708 ! Andrei's new change of 4.08.11                                (end)
2709 !**** ***************************************************************
2710
2711       subroutine output(sum_zh,sum_zv,sum_ldr,sum_kdp,sum_rhv,&
2712                                         sum_cdr,sum_ah,sum_adp,out)
2713
2714 !use microprm
2715
2716       !implicit double precision (a-h,o-z)
2717           implicit none
2718
2719           intrinsic dimag, datan2, dreal
2720 ! ### Interface
2721       double precision :: out(10), sum_zh, sum_zv, sum_ldr, sum_kdp, &
2722                                                   sum_cdr, sum_ah, sum_adp
2723       complex(8) :: sum_rhv
2724       double precision,parameter :: pi = 3.14159265D0
2725 ! ### Interface
2726
2727 !**** ************
2728 ! ZH output &
2729 !**** ************
2730
2731       if(sum_zh.gt.1.0d-2) then
2732         out(1) = 1.0d1*dlog10(sum_zh)
2733       else
2734         out(1) = -35.0d0
2735       endif
2736
2737 !**** ************
2738 ! ZV output &
2739 !**** ************
2740
2741       if(sum_zh.gt.1.0d-2) then
2742         out(2) = 1.0d1*dlog10(sum_zv)
2743       else
2744         out(2) = -10.0d0
2745       endif
2746
2747 !**** ************
2748 ! ZDR output &
2749 !**** ************
2750
2751       if(sum_zh.gt.1.0d-2) then
2752         out(3) = 1.0d1*dlog10(sum_zh/sum_zv)
2753       else
2754 ! Andrei's new change of 27.07.11                             (start)
2755         out(3) = -10.0d0
2756 !        out(3) = 0.0d0
2757 ! Andrei's new change of 27.07.11                               (end)
2758       endif
2759
2760 !**** ************
2761 ! LDR output &
2762 !**** ************
2763
2764       if(sum_zh.gt.1.0d-2) then
2765         !test_sum_zh = abs(sum_zh)
2766         if ( sum_zh.lt.(sum_ldr*10.0D10) ) then
2767           out(4) = 1.0d1*dlog10(sum_ldr/sum_zh)
2768         else
2769           out(4) = -99.9d9
2770         endif
2771       else
2772         out(4) = -100.0d0
2773       endif
2774
2775 !**** ************
2776 ! KDP output &
2777 !**** ************
2778
2779       if(sum_zh.gt.1.0d-2) then
2780 !        out(5) = max(-100.0D0,min(100.0D0,sum_kdp))
2781         out(5) = sum_kdp
2782       else
2783        out(5) = 0.0D0
2784       endif
2785
2786
2787 !**** ************
2788 ! RHV output &
2789 !**** ************
2790
2791       if(sum_zh.gt.1.0d-2) then
2792         out(6) = abs(sum_rhv)/sqrt(sum_zh*sum_zv)
2793       else
2794         out(6) = 0.0d0
2795       endif
2796
2797 !**** ***********
2798 ! DELTA output &
2799 !**** ***********
2800
2801       if(sum_zh.gt.1.0d-2) then
2802           out(7) = datan2(dimag(sum_rhv),dreal(sum_rhv))*180.0d0/pi
2803       else
2804         out(7) = 0.0d0
2805       endif
2806
2807 !**** ***********
2808 ! Simulated CDR output &
2809 !**** ***********
2810
2811       if(sum_cdr.gt.1.0d-8) then
2812         out(8) = 1.0d1*dlog10(sum_cdr)
2813       else
2814         out(8) = 0.0d0
2815       endif
2816
2817 !**** ***********
2818 ! AH output &
2819 !**** ***********
2820
2821       if(sum_zh.gt.1.0d-2) then
2822         out(9) = sum_ah
2823       else
2824         out(9) = 0.0d0
2825       endif
2826
2827
2828 !**** ***********
2829 ! ADP output &
2830 !**** ***********
2831
2832       if(sum_zh.gt.1.0d-2) then
2833         out(10) = sum_adp
2834       else
2835         out(10) = 0.0d0
2836       endif
2837
2838
2839       return
2840       end subroutine output
2841 ! subroutine output &
2842 ! ********************************************************************
2843           subroutine calc_scattering_snow (wl,snow_mass,den_dry,den_wet, &
2844                                                                    fract_mass_water,dc_water,dc_snow,dc_wet,f_a,f_b,f_a0,f_b0)
2845
2846       USE scatt_tables,ONLY:twolayer_snow,rpquada,usequad
2847       !USE t_matrix2_quad_mod                     ! ### [KS] : this is not linked since we use look-up-tables
2848       !USE t_matrix2_double_mod                   ! ### [KS] : this is not linked since we use look-up-tables
2849
2850       implicit none
2851
2852           intrinsic DCONJG
2853           ! ### Interface
2854           double precision :: wl, snow_mass, den_dry, den_wet, fract_mass_water
2855           complex(8) :: dc_water,dc_snow,dc_wet,f_a, f_b, f_a0, f_b0
2856           ! ### Interface
2857
2858           ! ### Local
2859           double precision :: pi, den_water
2860       double precision :: degree, dd_dry, dd_melt, dd, fvw, aspect_dry, aspect_melt, aspect
2861       double precision :: rp, rpquad, aspect2, dcore, angle, dd2, ff, ff2, shape2_a, shape2_b
2862       double precision :: shape1_a, shape1_b, psi, dd1, tmp, shape_a, shape_b
2863       complex(8) :: num, denum
2864       integer, parameter :: twolayer_snow_rayleigh = 1  ! ### [KS] : Always equal to 1 ?
2865           ! ### Local
2866
2867       data pi, den_water /3.14159265D0, 1.0d0/
2868
2869       degree = 1.0d0/3.0d0
2870
2871       dd_dry = 1.d1*(snow_mass/den_dry)**degree
2872
2873       dd = 1.d1*((6.0D0/pi)*snow_mass*(fract_mass_water/den_water+ &
2874            (1.0d0-fract_mass_water)/den_dry))**degree
2875
2876       fvw = den_dry*fract_mass_water/((1-fract_mass_water)*den_water+fract_mass_water*den_dry)
2877
2878       if(dd_dry.lt.1.0d1) then
2879          aspect_dry = 1.0d0-2.0d-2*dd_dry
2880       else
2881          aspect_dry = 0.8d0
2882       end if
2883
2884       dd_melt = 1.d1*((6.0D0/pi)*snow_mass/den_water)**degree
2885
2886       if(dd_melt.lt.1.0d1) then
2887          aspect_melt=0.9951d0+0.0251d0*dd_melt-0.03644*dd_melt**2+ &
2888                      0.005303*dd_melt**3-0.0002492*dd_melt**4
2889       else
2890          aspect_melt  = 0.4131d0
2891       end if
2892
2893       aspect = aspect_dry+fract_mass_water*(aspect_melt-aspect_dry)
2894
2895       if(aspect.eq.1.0d0) aspect=0.9999d0
2896
2897       rp = dd*abs(sqrt(dc_wet))/(1.d1*wl)
2898
2899       if(rp.gt.0.1d0) then
2900          if (twolayer_snow == 1) then
2901              aspect2 = aspect
2902              dcore=(1.d0-fvw)**(1.d0/3.d0)*dd
2903              rpquad = rpquada(3)
2904              if ((rp.lt.rpquad) .OR. (usequad .EQV. .FALSE.)) then
2905                  !call t_matrix2_dp(wl,dd,dcore,aspect,aspect2,dc_water,dc_snow,f_a,f_b,f_a0,f_b0)
2906                                          ! [KS] >> look-up-tables, this part is not reached
2907              else
2908                  !call t_matrix2_qp(wl,dd,dcore,aspect,aspect2,dc_water,dc_snow,f_a,f_b,f_a0,f_b0)
2909                                          ! [KS] >> look-up-tables, this part is not reached
2910              endif
2911
2912          else
2913              angle = 1.8d2
2914              ! call t_matrix(dd,wl,dc_wet,aspect,angle,f_a,f_b,6,'Snow') ! [KS] >> This is not linked as we use lookup tables
2915              angle = 0.0d0
2916              ! call t_matrix(dd,wl,dc_wet,aspect,angle,f_a0,f_b0,6,'Snow') ! [KS] >> This is not linked as we use lookup tables
2917              f_b0 = -DCONJG(f_b0)
2918              f_a0 = DCONJG(f_a0)
2919          endif
2920       else
2921          if (twolayer_snow_rayleigh == 1) then
2922 !     **************   external spheroid    ***************
2923
2924             dd2=1.d1*((6.0D0/pi)*snow_mass*(fract_mass_water/den_water+(1-fract_mass_water)/den_dry))**degree
2925             ff       = sqrt((1.0d0/aspect)**2-1.0d0)
2926             ff2      = ff**2
2927             shape2_a = ((1.0+ff2)/ff2)*(1.0d0-atan(ff)/ff)
2928             shape2_b = 0.5d0*(1.0d0-shape2_a)
2929
2930 !     ************   inner spheroid   ********************
2931             dd1=1.d1*((6.0D0/pi)*snow_mass*(1.0d0-fract_mass_water)/den_dry)**degree
2932             shape1_a = shape2_a
2933             shape1_b = shape2_b
2934             psi      = (dd1/dd2)**3
2935
2936 !     ******   Scattering amplitudes    *********
2937             tmp   =  pi**2*dd2**3/(6.0d2*wl**2)
2938             num   = (dc_water-1.0d0)*(dc_water+(dc_snow-dc_water)* &
2939                     (shape1_a-psi*shape2_a))+ &
2940                     psi*dc_water*(dc_snow-dc_water)
2941
2942             denum = (dc_water+(dc_snow-dc_water)*(shape1_a-psi*shape2_a))* &
2943                     (1.0d0+(dc_water-1.0d0)*shape2_a)+ &
2944                     psi*shape2_a*dc_water*(dc_snow-dc_water)
2945
2946             f_a0   = tmp*num/denum
2947             f_a  = DCONJG(f_a0)
2948
2949             num   = (dc_water-1.0d0)*(dc_water+(dc_snow-dc_water)* &
2950                     (shape1_b-psi*shape2_b))+ &
2951                      psi*dc_water*(dc_snow-dc_water)
2952
2953             denum=(dc_water+(dc_snow-dc_water)*(shape1_b-psi*shape2_b))* &
2954                   (1.0d0+(dc_water-1.0d0)*shape2_b)+ &
2955                    psi*shape2_b*dc_water*(dc_snow-dc_water)
2956             f_b0  = tmp*num/denum
2957             f_b = DCONJG(f_b0)
2958          else
2959             ff      = sqrt((1.0d0/aspect)**2-1.0d0)
2960             ff2     = ff**2
2961
2962             shape_a = ((1.0+ff2)/ff2)*(1.0d0-atan(ff)/ff)
2963
2964             shape_b = 0.5d0*(1.0d0-shape_a)
2965
2966             tmp     = pi**2*dd**3/(6.0d2*wl**2)
2967
2968             f_a0     = tmp/(shape_a+1.0d0/(dc_wet-1.0d0))
2969             f_b0     = tmp/(shape_b+1.0d0/(dc_wet-1.0d0))
2970             f_a    = dconjg(f_a0)
2971             f_b    = dconjg(f_b0)
2972          endif
2973       endif
2974
2975       return
2976       end subroutine calc_scattering_snow
2977 ! subroutine calc_scattering_snow &
2978 ! ****************************************
2979
2980 END MODULE module_mp_SBM_polar_radar
2981 !**** **************************************************************** &
2982 #endif