!***************************************************
! AGED BIOMASS-BURNING DPM AND SIGMA PARAMETERIZATION
!****************************************************

! Written by Kim Sakamoto, November 2015

! PURPOSE: The aged Dpm and sigma parameterization 
! uses an emulator of the SAM-TOMAS plume model to predict
! the coagulational growth of a single-mode biomass-burning 
! aerosol size distribution. 

! The 7 inputs to the parameterization are described (with ranges)
! below in the INPUTS section. This code returns estimated Dpm and
! sigma of the aged lognormal size distribution.  

! Further information can be found in:
! K. M. Sakamoto, R. G. Stevens, and J. R. Pierce,
! The evolution of biomass-burning aerosol size distributions 
! due to coagulation: dependence on fire and meteorological 
! details and parameterization., 
! Geosci. Model Dev. Discuss., 2015 

!--------- USAGE--------------------!

! call READ_FILES to load emulator matrix files 
! from dpm_emulator_files and sig_emulator_files directories
! must be called BEFORE param_dpm_sig

! call param_dpm_sig with 7 input parameters and two output 
! parameters

! EXAMPLE: 

! call READ_FILES()
! call param_dpm_sig(idpmemis=dpmemis, imassemis=massemis, isigemis=sigemis,&
!                    ifire_area=fire_area, iwindspeed=windspeed, itime=time, idmix=dmix &
!                    odpmout=dpmout, osigmaout=sigmaout)

!------------------------------------!
module param_BB

implicit none
real, allocatable :: xmod(:,:), ainv(:,:), g(:), ainvh(:,:), inv_hainvh(:,:)
real, allocatable :: h(:), betahat(:), roughness(:), min_d(:), max_d(:), t(:), ainvt(:)
real, allocatable :: tainvh(:), hdiff(:), xpred(:)

real, allocatable :: xmod_s(:,:), ainv_s(:,:), g_s(:), ainvh_s(:,:), inv_hainvh_s(:,:)
real, allocatable :: h_s(:), betahat_s(:), roughness_s(:), min_s(:), max_s(:), t_s(:), ainvt_s(:)
real, allocatable :: tainvh_s(:), hdiff_s(:), xpred_s(:)

real prec, scalemean, scalesd
integer nmod, nreg, ninputs

integer nmod_s, nreg_s, ninputs_s
real prec_s, scalemean_s, scalesd_s

contains

  SUBROUTINE param_dpm_sig(idpmemis, imassemis, isigemis, ifire_area, &
                           iwindspeed, itime, idmix, &
                           odpmout, osigout)
  ! Returns dpm and sigma of the aged BB size distribution based on fire
  ! and meteorological conditions

    !------INPUTS------------------------------------------------------
    
    real, intent(in):: idpmemis
    ! initial emitted median diameter [nm]
    ! tested range 20-100 nm
    
    real, intent(in):: imassemis
    ! mass emissions flux [mg m-2 s-1]
    ! tested range 0.02-5 mg m-2 s-1

    real, intent(in):: isigemis
    ! initial modal width of distribution [unitless]
    ! tested range 1.2-2.4

    real, intent(in):: ifire_area
    ! fire area [km2]
    ! tested range 1-49 km2

    real, intent(in):: iwindspeed
    ! boundary layer windspeed [ms-1]
    ! tested range 2-20 ms-1

    real, intent(in):: itime
    ! time since emission [min]
    ! tested range 0-300 min

    real, intent(in):: idmix
    ! mixing depth of aerosol layer [m]
    ! tested range: 120-2500 m

    !-----OUTPUTS-------------------------------------------------------
    
    real, intent(out):: odpmout
    ! aged median diameter of lognormal size distribution [nm]

    real, intent(out):: osigout
    ! aged modal width of lognormal size distribution [unitless]

    !-----Reading inputs-------------------
   
    real :: dpmemis, massemis, sigemis, fire_area, windspeed, time, dmix
    real sigout, dpmout

    dpmemis=idpmemis
    massemis=imassemis
    sigemis=isigemis
    fire_area=ifire_area
    windspeed=iwindspeed
    time=itime
    dmix=idmix
    
    !-----Code--------------------------------
    
    odpmout=dpm_emulator(dpmemis, massemis,sigemis,fire_area,windspeed,time, dmix)
    osigout=sig_emulator(dpmemis,massemis,sigemis,fire_area, windspeed, time, dmix)

    END SUBROUTINE param_dpm_sig

!-----FUNCTION DEFINITIONS----------------------------------------------------------
    real function sig_emulator(dpmemis, massemis, sigemis, fire_area, windspeed, time, dmix)
    ! evaluate the dpm emulator at the parameter point
    !------------INPUTS----------------------
    real :: dpmemis, massemis, sigemis, fire_area, windspeed, time, dmix
    integer i

    !-------------CODE-----------------

    ! put certain parameters in 1og10 format as required (dpm, massemiss, area and mixing depth)
    xpred_s(1)=log10(dpmemis)
    xpred_s(2)=log10(massemis)
    xpred_s(3)=sigemis
    xpred_s(4)=log10(fire_area)
    xpred_s(5)=windspeed
    xpred_s(6)=time
    xpred_s(7)=log10(dmix)

    ! rescale inputs exactly as in the emulator building code
    xpred_s = (xpred_s - min_s)/(max_s-min_s)

    h_s(1) = 1.
    IF (nreg_s.gt.1) THEN
    DO i=1, ninputs_s
    h_s(i+1) = xpred_s(i)  ! assumes linear regression terms only
    ENDDO
    ENDIF

    ! compute t vector of correlations between xpred and training data inputs
    DO i=1, nmod_s
      t_s(i) = exp(-SUM(roughness_s*(xpred_s-xmod_s(i,:))*(xpred_s-xmod_s(i,:))))
    ENDDO

    ! compute mean
    sig_emulator = scalemean_s + scalesd_s*(DOT_PRODUCT(h_s, betahat_s) + DOT_PRODUCT(t_s, g_s))

    return
    end function sig_emulator

!------------------------------------------------------------
    real function dpm_emulator(dpmemis, massemis, sigemis, fire_area, windspeed, time, dmix)
    ! evaluate the dpm emulator at the parameter point
    !use emulator
    !------------INPUTS----------------------
    real :: dpmemis, massemis, sigemis, fire_area, windspeed, time, dmix
    integer i
    !-------------CODE-----------------

    ! put certain parameters in 1og10 format as required (dpm, massemiss, area and mixing depth)
    xpred(1)=log10(dpmemis)
    xpred(2)=log10(massemis)
    xpred(3)=sigemis
    xpred(4)=log10(fire_area)
    xpred(5)=windspeed
    xpred(6)=time
    xpred(7)=log10(dmix)

    ! rescale inputs exactly as in the emulator building code
    xpred = (xpred - min_d)/(max_d-min_d)

    h(1) = 1.
    IF (nreg.gt.1) THEN
    DO i=1, ninputs
      h(i+1) = xpred(i)  ! assumes linear regression terms only
    ENDDO
    ENDIF

    ! compute t vector of correlations between xpred and training data inputs
    DO i=1, nmod
      t(i) = exp(-SUM(roughness*(xpred-xmod(i,:))*(xpred-xmod(i,:))))
    ENDDO

    ! compute mean
    dpm_emulator = scalemean + scalesd*(DOT_PRODUCT(h, betahat) + DOT_PRODUCT(t, g))
    dpm_emulator=10**dpm_emulator ! emulated dpm in [nm]

    return
    end function dpm_emulator

!------------------------------------------------------------
! READ IN EMULATOR FILES
!-----------------------------------------------------------

    SUBROUTINE READ_FILES()
    ! open emulator files (both for dpm and sigma),     read data
    implicit none
    integer i, j

    !open dpm emulator setup files
    OPEN(UNIT=11, FILE='dpm_emulator_files/emulator_ainv_dpm.txt', status='old')  ! inverse correlation matrix
    OPEN(UNIT=12, FILE='dpm_emulator_files/emulator_training_inputs_dpm.txt', status='old') ! training data inputs
    OPEN(UNIT=13, FILE='dpm_emulator_files/emulator_mu_out_dpm.dat', status='old') ! estimated regression parameters
    OPEN(UNIT=14, FILE='dpm_emulator_files/emulator_precision_out_dpm.dat', status='old') ! estimated GP prec = 1/variance
    OPEN(UNIT=15, FILE='dpm_emulator_files/emulator_ainvh_dpm.txt', status='old') ! ainv.H
    OPEN(UNIT=16, FILE='dpm_emulator_files/emulator_rough_out_dpm.dat', status='old') ! estimated function roughnesses
    OPEN(UNIT=17, FILE='dpm_emulator_files/emulator_g_dpm.txt', status='old') ! ainv.(y-Hb)
    OPEN(UNIT=18, FILE='dpm_emulator_files/emulator_inv_hainvh_dpm.txt', status='old') ! (H'.ainv.H)^{-1}
    OPEN(UNIT=19, FILE='dpm_emulator_files/emulator_minmax_dpm.txt', status='old') ! maxs and mins of each input
    OPEN(UNIT=20, FILE='dpm_emulator_files/emulator_scale_dpm.txt', status='old') ! mean and sd of output scaling

    read(12,*)nmod, ninputs, nreg
    ALLOCATE(h(nreg), t(nmod), min_d(ninputs), max_d(ninputs), ainv(nmod,nmod), xmod(nmod,ninputs), &
     ainvh(nmod,nreg), betahat(nreg), roughness(ninputs), g(nmod), inv_hainvh(nreg,nreg), &
     ainvt(nmod),tainvh(nreg), hdiff(nreg), xpred(ninputs))

    read(19,*)(min_d(i), i=1, ninputs)
    read(19,*)(max_d(i), i=1, ninputs)
    read(20,*)scalemean, scalesd
    DO i = 1, nmod
      read(11,*)(ainv(i,j), j = 1, nmod)
      read(12,*)(xmod(i,j), j = 1, ninputs)
      read(15,*)(ainvh(i,j), j = 1, nreg)
    ENDDO

    read(13,*)(betahat(i), i=1, nreg)
    read(14,*)prec
    read(16,*)(roughness(i), i=1, ninputs)
    read(17,*)(g(i), i=1, nmod)

    DO i = 1, nreg
      read(18,*)(inv_hainvh(i, j), j = 1, nreg)
    ENDDO

    DO i=11,20
       CLOSE(i)
    ENDDO

! SIGMA EMULATOR
    OPEN(UNIT=21, FILE='sig_emulator_files/emulator_ainv_sig.txt', status='old')  ! inverse correlation matrix
    OPEN(UNIT=22, FILE='sig_emulator_files/emulator_training_inputs_sig.txt', status='old') ! training data inputs
    OPEN(UNIT=23, FILE='sig_emulator_files/emulator_mu_out_sig.dat', status='old') ! estimated regression parameters
    OPEN(UNIT=24, FILE='sig_emulator_files/emulator_precision_out_sig.dat', status='old') ! estimated GP prec = 1/variance
    OPEN(UNIT=25, FILE='sig_emulator_files/emulator_ainvh_sig.txt', status='old') ! ainv.H
    OPEN(UNIT=26, FILE='sig_emulator_files/emulator_rough_out_sig.dat', status='old') ! estimated function roughnesses
    OPEN(UNIT=27, FILE='sig_emulator_files/emulator_g_sig.txt', status='old') ! ainv.(y-Hb)
    OPEN(UNIT=28, FILE='sig_emulator_files/emulator_inv_hainvh_sig.txt', status='old') ! (H'.ainv.H)^{-1}
    OPEN(UNIT=29, FILE='sig_emulator_files/emulator_minmax_sig.txt', status='old') ! maxs and mins of each input
    OPEN(UNIT=30, FILE='sig_emulator_files/emulator_scale_sig.txt', status='old') ! mean and sd of output scaling

    read(22,*)nmod_s, ninputs_s, nreg_s

    ALLOCATE(h_s(nreg_s), t_s(nmod_s), min_s(ninputs_s),&
    max_s(ninputs_s), ainv_s(nmod_s,nmod_s),&
    xmod_s(nmod_s,ninputs_s), ainvh_s(nmod_s,nreg_s),&
    betahat_s(nreg_s), roughness_s(ninputs_s),&
    g_s(nmod_s), inv_hainvh_s(nreg_s,nreg_s),&
    ainvt_s(nmod_s),tainvh_s(nreg_s),&
    hdiff_s(nreg_s), xpred_s(ninputs_s))


    read(29,*)(min_s(i), i=1, ninputs_s)
    read(29,*)(max_s(i), i=1, ninputs_s)
    read(30,*)scalemean_s, scalesd_s

    DO i = 1, nmod_s
      read(21,*)(ainv_s(i,j), j = 1, nmod_s)
      read(22,*)(xmod_s(i,j), j = 1, ninputs_s)
      read(25,*)(ainvh_s(i,j), j = 1, nreg_s)
    ENDDO

    read(23,*)(betahat_s(i), i=1, nreg_s)
    read(24,*)prec_s
    read(26,*)(roughness_s(i), i=1, ninputs_s)
    read(27,*)(g_s(i), i=1, nmod_s)

    DO i = 1, nreg_s
      read(28,*)(inv_hainvh_s(i, j), j = 1, nreg_s)
    ENDDO

    DO i=20,30
      CLOSE(i)
    ENDDO

    END SUBROUTINE read_files

end module param_BB
!-----------END CODE ----------------------------!
