program predict_psd


c This procedure uses the second moment and in-cloud temperature
c to predict the ice particle size distribution. 
c The moment relations use the 2nd moment as a reference moment.
c If IWC is computed using a different moment (y ne 2) then
c the second moment needs to be estimated first.
c (To avoid this step the moment relations could be found for a different
c reference moment - requires reanalysis of the original data).

c For more details please read "Parametrization of ice particle size
c distributions for mid-latitude stratiform cloud" by Field et al. 2005 QJRMS 1997-2017 

c !!!!!!!!!!!variables!!!!!!!!!!!!!!!!!!!!
c Input
c tc in C
c iwc in g m^-3
c x,y prefactor and exponent in mass-size relation m=xD^y
c for the UM    x=0.069, y=2.0    (SI units)

c Output
c D particle size in microns
c dN_dD particle size distribution: dN/dD m^-4
c L23 ratio of M3/M2 - mean mass weighted size if y=2


       real D(500),dN_dD(500),phi(500),a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10
       real x,y,tc,iwc,iwcsi,My,m2,n,loga_,a_,b_,m3,xx


c!!!constants for moment relations based on 2nd moment as reference!!!!!!!!!!!
       a1=5.065339
       a2=-0.062659
       a3=-3.032362
       a4=0.029469
       a5=-0.000285
       a6=0.31255
       a7=0.000204
       a8=0.003199
       a9=0.0
       a10=-0.015952

       b1=0.476221
       b2=-0.015896
       b3=0.165977
       b4=0.007468
       b5=-0.000141
       b6=0.060366
       b7=0.000079
       b8=0.000594
       b9=0.0
       b10=-0.003577
c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
c    read in input file
       open(unit=2,file="input.txt")
       read(2,*) iwc,tc,x,y
       close(2)





c incloud temperature check

      if (tc.gt.0) then
         print*,'in-cloud temperature too warm (>0C)'
         goto 1000
      endif  
      if (tc.lt.-60) then 
         print*,'Warning: in-cloud temp colder than in original data'
      endif  

c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!



c Calculate 2nd moment
      iwcsi=iwc/1e3  !kg m^-3

      My=iwcsi/x  !yth moment

c use moment relations to get 2nd moment required for predicting psd

      m2=My 
c default for UM y=2 (do not need to do inversion in following if-block to avoid slight error introduced by plane fitting)


      if (y.ne.2) then !estimate 2nd moment if y ne 2
        n=y
        loga_=a1+a2*tc+a3*n+a4*tc*n+a5*tc**2+a6*n**2+a7*tc**2*n+
     $     a8*tc*n**2+a9*tc**3+a10*n**3
        a_=10**loga_
        b_=b1+b2*tc+b3*n+b4*tc*n+b5*tc**2+b6*n**2+b7*tc**2*n+
     $     b8*tc*n**2+b9*tc**3+b10*n**3
        M2=(My/a_)**(1.0/b_)
      endif


c !!M2 range check (1e-5 3e-2 m^-1)!!!!!!!!!!!!!!!!!!!!!

      if (m2.lt.0.0) then
         print*,'Negative M2'
         goto 1000
      endif
      if (m2.lt.1e-5.or.m2.gt.3e-2) then
        print*,'Warning: M2 outside of range in original data'
      endif 
c !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


c !Use 2nd and 3rd moment to predict psd!!!!!!!!!!!!!!!!!!!!


      n=3
        loga_=a1+a2*tc+a3*n+a4*tc*n+a5*tc**2+a6*n**2+a7*tc**2*n+
     $     a8*tc*n**2+a9*tc**3+a10*n**3
        a_=10**loga_
        b_=b1+b2*tc+b3*n+b4*tc*n+b5*tc**2+b6*n**2+b7*tc**2*n+
     $     b8*tc*n**2+b9*tc**3+b10*n**3
      m3=a_*m2**b_



c !!!define universal psd!!!!
      do 10 i=1,500
        xx=i/40. !dimensionless size
        phi(i)=490.6*exp(-20.78*xx)+17.46*xx**0.6357*exp(-3.29*xx) !universal distribution
c	print*,xx,phi(i),m2,m3,a_,b_
        D(i)=1e6*xx*(m3/m2)
        dN_dD(i)=phi(i)*m2**4/m3**3
 10   continue




c !write output file

      open(unit=2,file="output.txt")
      write(2,*) 'Input values'
      write(2,*) 'IWC [g/m3]=',iwc
      write(2,*) 'T [C] =',tc
      write(2,*) 'x , y [SI]=', x, y
      write(2,*) ' '
      write(2,*) 'Output values'
      write(2,*) 'L32 [microns]=',1e6*m3/m2
      write(2,*) 'D (microns) dN/dD (m^-4)'
      do 20 i=1,500
        write(2,*),D(i),dN_dD(i)
 20   continue
      close(2)	       
 

 1000   continue

        stop 
        end