ss_model.f90

SUBROUTINE SS_MODEL(CTRL,FREQ,TETAD,Y,TB_UBC,AUX_INS,TB_OUT,pixel,replicate,&
  rank,meas) 

! -----------------------------------------------------------------------------
!
!     SOIL AND SNOW RADIATIVE TRANSFER MODEL
!
!     FORTRAN 90 VERSION OF MEMLS PROGRAM (ORIGINALLY WRITTEN IN MATLAB BY 
!     WIESMANN AND MATZLER TO COMPUTE BRIGHTNESS TEMPERATURES OF N-LAYERED 
!     SNOWPACKS) INTEGRATED WITH A SOIL RADIATIVE TRANSFER SCHEME.  THIS 
!     VERSION OF THE CODE IS SET UP TO PERFORM CALCULATIONS BASED ON DATA FROM
!     A SINGLE SNOWPACK AT THE FREQUENCIES AND ANGLES SPECIFIED IN THE INPUTS.  
!     IF THE NUMBER OF LAYERS SPECIFIED IN CTRL(2) IS ZERO, THEN THIS CODE
!     WILL COMPUTE THE BARE SOIL BRIGHTNESS TEMPERATURE
!
!     IN THE CURRENT VERSION, EACH SNOWPACK MUST HAVE THE SAME NUMBER OF LAYERS.
!     IN ADDITION, THIS CODE IS DESIGNED TO WORK WITH SNOWPACKS WITH INCOHERENT
!     LAYERS OR WITH ONE COHERENT LAYER.  MORE SPECIFICALLY, THE CODE
!     CANNOT HANDLE THE CASE WHERE THERE ARE SUCCEEDING COHERENT LAYERS.  IN 
!     THIS CASE, THE PROGRAM WILL WRITE 999 TO THE OUTPUT FILE.  
!
! ----------------------------------------------------------------------------
!
!     DESCRIPTION OF INPUTS AND OUTPUTS (SEE JOURNAL ENTRY 19 SEP 05)
!
!     THE CTRL VECTOR SPECIFIES NINE OVERALL RUN CONTROLS, NOT ALL OF WHICH 
!     ARE USED IN THIS MODEL.  CTRL=(N_CACLS,N_LYRS,ATM_SWITCH,CAN_SWITCH,
!       N_AUX_INS,N_SNOW_INS,N_CAN_INS,N_ATM_INS,N_FREQ)
!     THE FREQ AND TETAD VECTORS SPECIFY AT WHICH FREQUENCIES CALCULATIONS 
!       ARE TO BE DONE.
!     THE Y ARRAY SPECIFIES SNOW INPUTS.  THERE IS ONE ROW FOR EACH SNOW
!       LAYER SPECIFIED IN CTRL(2), ONE COLUMN SPECIFIED FOR EACH SNOW
!       VARIABLE,IN ORDER: LAYER THICKNESS [M], LAYER DENSITY [KG/M3],
!       LAYER GRAIN DIAMETER [M], LAYER LIQUID WATER CONTENT [FRAC],
!       LAYER TEMPERATURE [K].
!     THE TB_UBC ARRAY SPECIFIES THE UPPER BOUNDARY CONDITION BRIGHTNESS
!       TEMPERATURE (SKY, VEGETATION, COSMIC).  IT HAS TWO ROWS (FOR H AND V
!       POLARIZATION) AND CTRL(9) COLUMNS: ONE FOR EACH FREQUENCY.
!     THE AUX_INS ARRAY SPECIFIES SOME AUXILIARY INPUTS.  THESE INCLUDE (IN 
!       ORDER): N_LYRS, SOIL TEMPERATURE [K], SOIL SATURATION [FRAC], SOIL
!       POROSITY [FRAC], AND THE CONSTANT OF PROPORTIONALITY BETWEEN THE 
!       GRAIN SIZE AND THE CORRELATION LENGTH.
!
!     THE TB_OUT ARRAY CONTAINS THE CALCULATED BRIGHTNESS TEMPERATURE OUTPUTS
!       AT EACH POLARIZATION (ROWS) AND FREQUENCY (COLUMNS).
!            
! ----------------------------------------------------------------------------
!
!     1.DESCRIPTION.  THE CODE IS A TRANSLATION OF THE MEMLS2 CODE RECEIVED 
!     THROUGH CORRESPONDENCE WITH MATZLER AND A SOIL RADIATIVE TRANSFER
!     MODEL RECEIVED THROUGH CORRESPONDENCE WITH PULLIAINEN SOMETIME IN 2003.
!     THE MEMLS CODE USES THE GROUND REFLECTIVITY AND THE SNOWPACK DATA TO 
!     COMPUTE THE BRIGHTNESS TEMPERATURE OF THE SNOWPACK. THE CODE WAS 
!     TRANSLATED TO FORTRAN IN ORDER  TO SAVE COMPUTATION TIME IN A DATA 
!     ASSIMILATION SCHEME IN WHICH MILLIONS OF CALLS MUST BE MADE TO THE 
!     RADIATIVE TRANSFER MODEL .  IN TESTS, ONE HUNDRED BRIGHTNESS TEMPERATURE
!     CALCULATIONS TOOK TEN TIMES LESS PROCESSOR TIME IN THIS FORTRAN VERSION 
!     THAN IN THE ORIGINAL MATLAB.  EXTENSIVE VERIFICATION WAS PERFORMED USING
!     SNOWPIT DATA FOR NINE DAYS IN ONE WINTER DURING THE CLPX PROGRAM AS WELL
!     AS 3-LAYER SNOWMODEL RESULTS FROM THE MODIFIED SAST (SUN ET AL, 99 IN 
!     IN JGR) + SSIB MODEL RUN AT MAMMOTH MOUNTAIN, TO ENSURE THAT THE FORTRAN
!     CODE PRODUCES IDENTICAL RESULTS AS THE MATLAB CODE.
!
!     2. MODIFICATION.  THE MAIN CHANGE MADE TO THE MEMLS CODE IN THE 
!     TRANSLATION PROCESS WAS THAT THE VAN POLDER APPROXIMATION (EFFECTIVE 
!     MIXING THEORY) WAS  NO LONGER BEING USED IN THE MATLAB CODES WHEN I 
!     RECEIVED THEM.  IN FACT, THE REAL PART OF THE SNOW DIELECTRIC WAS USED 
!     INSTEAD IN SEVERAL SUBROUTINES. THIS WAS SLIGHTLY INCONSISTENT WITH THE 
!     PAPERS, SO (FOLLOWING MATZLER'S ADVICE VIA EMAIL) I WROTE A NEW 'POLDER' 
!     SUBROUTINE WHICH USES THE NEWTON-RAPHSON APPROXIMATION TO SOLVE FOR THE 
!     SNOW EFFECTIVE PERMITTIVITY - SEE MATZLER'S 1996 PAPER IN IEEE FOR THE 
!     EQUATIONS AND THE SUBROUTINE COMMENTS BELOW.  ACCORDING TO COMPARISON 
!     WITH THE SNOWPIT DATA RESULTS, THE EFFECT OF USING THE SNOW DIELECTRIC 
!     INSTEAD OF EFFECTIVE PERMITTIVITY ON THE BRIGHTNESS TEMPREATURE WAS 
!     MINIMAL (~0.1 DEGREE).  HOWEVER, BECAUSE I COULD NOT TEST WITH A MORE 
!     EXTENSIVE DATASET, I LEFT THE SLIGHTLY MORE COMPUTATIONALLY EXPENSIVE 
!     EFFECTIVE PERMITTIVITY CALCULATION IN PLACE FOR THE SAKE OF INTERNAL 
!     CONSISTENCY.
!
!     ANOTHER CHANGE THAT I MADE TO MEMLS WAS TO SET UP THE PROGRAM TO PASS IN 
!     SNOW GRAIN DIAMETER SINCE THIS IS USUALLY WHAT IS MEASUSRED IN SNOWPITS 
!     AND MODELED BY PROGNOSTIC EQUATIONS IN SNOW MODELS.  THE CONSTANT OF 
!     PROPORTIONALITY BETWEEN THE GRAIN DIAMETER AND THE CORRELATION LENGTH (SEE
!     MATLZER, 2002 IN JOURNAL OF GLACIOLOGY) IS ALSO PASSED IN.  
!
!     A THIRD CHANGE I MADE TO MEMLS WAS TO USE A SWITCH TO DETERMINE WHETHER 
!     THE BORN APPROXIMATION OR THE EMPIRICAL SCATTERING COEFFICIENT WOULD BE 
!     USED.  IF THE MAXIMUM CORRELATION LENGTH IN THE SNOWPACK IS GREATER THAN 
!     0.33 MM, THE BORN APPROXIMATION IS USED; OTHERWISE, THE EMPIRICAL FORMULA
!     IS USED.
!
!     3. COMPLETENESS.  I TRANSLATED EVERYTHING PRECISELY AS IT WAS IN THE 
!     MATLAB CODE AS CLOSELY AS POSSIBLE.  THE ONLY EXCEPTION WAS THAT (BECAUSE
!     I TYPICALLY DON'T WORK WITH SNOWPIT DATA WITH VERY THIN LAYERS) I REMOVED
!     THE ABILITY TO DEAL WITH SUCCEEDING COHERENT SNOWPACK LAYERS IN THE 
!     INTEREST OF FINISHING THE TRANSLATION IN A TIMELY MANNER.  I TRANSLATED 
!     ONLY THE BEGINNING OF THE SLRED SUBROUTINE, AND SET IT UP TO WRITE 999 TO
!     THE OUTPUT FILE IF IT FINDS SUCCEEDING COHERENT LAYERS.  ISOLATED COHERENT
!     LAYERS ARE HANDLED CORRECTLY. THE EXCEPTION IS FOR THREE-LAYER SNOWPACKS;
!     ANY COMBINATION OF COHERENT LAYERS SHOULD GIVE CORRECT RESULTS FOR THREE-
!     LAYER PACKS.
!     

IMPLICIT NONE

! 0.  PRELIMINARIES
!  A. DECLARATIONS
!   0.INPUTS AND OUTPUTS
INTEGER,INTENT(IN) :: CTRL(9)
REAL,INTENT(IN) :: AUX_INS(5),TB_UBC(2,CTRL(9)),FREQ(CTRL(9)),TETAD(CTRL(9))
REAL,INTENT(IN) :: Y(CTRL(2),CTRL(6))
REAL,INTENT(OUT):: TB_OUT(2,CTRL(9))

!   1.FOR PART 1
REAL, DIMENSION(:), ALLOCATABLE :: SNOWDATA,TI,WIFR,ROIKG,DI,GDI,&
ROI,GDIMM,PCI,EPSI,EPSII,GAI,NS,TEI,DEI,SIH,SIV
REAL,DIMENSION(:), ALLOCATABLE :: XNUM,XROI,XEPSI,XEPSII,XTEI,&
XSIH,XSIV,XDI,XDEI,XTI,XPCI,XWIFR,XGAI
REAL,DIMENSION(:), ALLOCATABLE :: RROI,REPSI,REPSII,RTEI,RSIH,&
RSIV,RDI,RDEI,RTI,RPCI,RWIFR,RGAI
REAL GND_SIG,GND_EPS,MJU0,EPS0,PI,TGS,SNGDPCI,TSKYH,&
TSKYV,SASTPOROS,TETA,SOILSATURATION,GND_TEMP,GND_MV,MAXPCI,&
VEG_HEIGHT,VEG_TEMP,MG,S,NDL_THK,NDL_LGTH,NDL_NUM,T0,P0,MOIST0
INTEGER SCCHO,I,NUM,NFREQ,K,RNUM,MONTH,VEG_SWITCH,ATM_SWITCH
!   2.FOR PART 2
REAL EPS_T,TETAD_SOIL,EPS_B
COMPLEX, DIMENSION(:), ALLOCATABLE :: EPSR
COMPLEX EPS_UPPER
REAL KSIG,S0H,S0V

!   3.FOR PART 3
REAL,DIMENSION(:),ALLOCATABLE :: GBIH,GBIV,GS6,GA2I,TSCAT,&
  RSIHLONG,RSIVLONG
!   4.FOR PART 4
REAL,DIMENSION(:),ALLOCATABLE :: RI,TRI,DH,DV
!REAL TSG_H,TSG_V,ESG_H,ESG_V

!   5.FOR PART 5
COMPLEX EVEG,ESW
REAL TCAN_V,TCAN_H,CAN_LOSS_H,CAN_LOSS_V,EFSG_H,EFSG_V,TFSG_H,&
  TFSG_V,TATM,ATM_TRAN,TAFSG_H,TAFSG_V

integer,intent(in):: pixel,replicate,rank,meas

!  B. CONTROL STATEMENTS 

NFREQ=CTRL(9)  ! THIS CODE LOOPS OVER THE NUMBER OF FREQUENCIES FOR ONE SET
               ! OF SNOW PROPERTIES.
!NUM=aux_ins(1) ! NUMBER OF LAYERS IN THE SNOWPACK 
NUM=ctrl(2) ! NUMBER OF LAYERS IN THE SNOWPACK 

!  C. ALLOCATE AND EXTRACT STATEMENTS

ALLOCATE(DI(1:NUM),ROIKG(1:NUM),GDI(1:NUM),WIFR(1:NUM),TI(1:NUM))

DI(1:NUM)=   Y(1:NUM,1)
ROIKG(1:NUM)=Y(1:NUM,2)
GDI(1:NUM)=  Y(1:NUM,3)
WIFR(1:NUM)= Y(1:NUM,4)
TI(1:NUM)=   Y(1:NUM,5)

GND_TEMP=AUX_INS(2)
SOILSATURATION=AUX_INS(3)
SASTPOROS=AUX_INS(4)
SNGDPCI=AUX_INS(5)

GND_MV=SOILSATURATION*SASTPOROS
MJU0=1.2566e-006
EPS0=8.8542e-012
!  1. CHECK TO MAKE SURE THERE IS SNOW.  OTHERWISE DEAL WITH SOIL ONLY 
!     NOTE: AUX_INS(1) CONTAINS 0 IN THE CASE OF NO SNOW, THOUGH NUM
!     WILL BE 1 TO AVOID INCORRECT ARRAY ALLOCATIONS
IF (AUX_INS(1).EQ.0) THEN
  DO K=1,NFREQ
    ! THERE IS NO SNOWPACK, SO COMPUTE REFLECTIVITY OF GROUND AND OBTAIN
    ! OUTPUT BRIGHTNESS TEMPERATURES
    EPS_T=1.
    TETAD_SOIL=TETAD(K)
    EPS_UPPER=1.0003 !DIELECTRIC PERMITTIVITY OF AIR
    gnd_sig=0.002 !i think...

    KSIG=REAL(2*PI*FREQ(K)*1E9*(MJU0*EPS0*EPS_UPPER)**0.5)*GND_SIG
    CALL RUFFSOIL(FREQ(K)*1E9,GND_MV,GND_TEMP,KSIG,TETAD_SOIL,&
      EPS_UPPER,S0H,S0V)

    TSKYH=TB_UBC(1,K)
    TSKYV=TB_UBC(2,K)

    TB_OUT(1,K)=(1-S0H)*GND_TEMP+S0H*TSKYH
    TB_OUT(2,K)=(1-S0V)*GND_TEMP+S0V*TSKYV
  END DO

  RETURN !TO MAIN PROGRAM

END IF

! 2. COMPUTE BRIGHTNESS TEMPERATURE OF SNOWPACK AND SOIL

DO K=1,NFREQ

  !  A. PRELIMINARIES 
  !  A.1. DEFINE CONSTANTS

  GND_SIG=0.002
  GND_EPS=0
!  MJU0=1.2566e-006
!  EPS0=8.8542e-012

  PI=3.14159


  !  A.2. ALLOCATE AND UNIT CONVERSION STATEMENTS
  ALLOCATE(ROI(1:NUM),GDIMM(1:NUM),PCI(1:NUM))
  
  ROI=ROIKG/1000
  TETA=TETAD(K)*PI/180
  GDIMM=GDI*1000
  PCI=GDIMM*SNGDPCI

  !  A.3. DETERMINE WHICH SCATTERING COEFFICIENT COMPUTATION METHOD TO 
  !  USE THE VALUE OF 0.33 USED HERE FROM PERSONAL CORRESPONDENCE WITH 
  !  MATZLER 

  MAXPCI=MAXVAL(PCI)
  IF (MAXPCI<0.33) THEN
    SCCHO=1 ! USE EMPIRICAL SCATTERING COEFFICIENT FROM W&M 99
  ELSE 
    SCCHO=2 ! USE BORN APPROXIMATION FROM MATZLER AND WIESMANN 99
  ENDIF

!THIS OPTION OUGHT TO BE AUTOMATED:
!DURING MY WORK WITH SSIB3+MEMLS TO MODEL GBMR-7 TB, I NEEDED TO USE SCCHO=2
SCCHO=2
!IN ORDER TO SAVE TIME FOR THE SYNTHETIC TESTS, HOWEVER, I USED SCCHO=1
!SCCHO=1

  ! B.  COMPUTE RADIATIVE TRANSFER PROPERTIES OF SNOW

  ALLOCATE(EPSI(1:NUM),EPSII(1:NUM),GAI(1:NUM),NS(1:NUM),&
    TEI(1:NUM+1),DEI(1:NUM),SIH(1:NUM),SIV(1:NUM),XROI(1:NUM),&
    XEPSI(1:NUM),XEPSII(1:NUM),XTEI(1:NUM+1),XSIH(1:NUM),&
    XSIV(1:NUM),XDI(1:NUM),XDEI(1:NUM),XTI(1:NUM),XPCI(1:NUM),&
    XWIFR(1:NUM),XGAI(1:NUM))

  CALL RO2EPSD(ROI,TI,FREQ(K),EPSI,EPSII,NUM,pixel,replicate,rank,meas,k)
  CALL MIXMOD(FREQ(K),TI,WIFR,EPSI,EPSII,NUM)
  CALL ABSCOEFF(EPSI,EPSII,TI,FREQ(K),WIFR,GAI,NUM)
  NS=EPSI**0.5
  TEI(1:NUM)=ASIN(SIN(TETA)/NS)
  TEI(NUM+1)=TETA
  CALL PFADI(TEI,DI,DEI,NUM)
  CALL FRESNELC(TEI,EPSI,SIH,SIV,NUM)

  CALL SLRED(NUM,ROI,EPSI,EPSII,TEI,SIH,SIV,DI,DEI,TI,PCI,WIFR,&
    GAI,FREQ(K),RNUM,XROI,XEPSI,XEPSII,XTEI,XSIH,XSIV,XDI,XDEI,&
    XTI,XPCI,XWIFR,XGAI,pixel,replicate,rank,meas)

  ALLOCATE(RROI(RNUM),REPSI(RNUM),REPSII(RNUM),&
    RTEI(RNUM+1),RSIH(RNUM),RSIV(RNUM),RDI(RNUM),&
    RDEI(RNUM),RTI(RNUM),RPCI(RNUM),RWIFR(RNUM),&
    RGAI(RNUM))

  DO I=1,RNUM
    RROI(I)=XROI(I)
    REPSI(I)=XEPSI(I)
    REPSII(I)=XEPSII(I)
    RSIH(I)=XSIH(I)
    RSIV(I)=XSIV(I)
    RDI(I)=XDI(I)
    RDEI(I)=XDEI(I)
    RTI(I)=XTI(I)
    RPCI(I)=XPCI(I)
    RWIFR(I)=XWIFR(I)
    RGAI(I)=XGAI(I)
  END DO
  DO I=1,RNUM+1
    RTEI(I)=XTEI(I)
  END DO

  ! C.  COMPUTE GROUND REFLECTIVITIES

  ALLOCATE(EPSR(1:RNUM))

  EPS_T=1.
  TETAD_SOIL=TETAD(K)

  DO I=RNUM,1,-1
    EPS_B=EPSI(I)
    CALL AOT(EPS_T,EPS_B,TETAD_SOIL)
    EPS_T=EPS_B
  END DO

  EPSR=CMPLX(REPSI,(-1*REPSII)) 
  EPS_UPPER=EPSR(1)

  KSIG=REAL(2*PI*FREQ(K)*1E9*(MJU0*EPS0*EPS_UPPER)**0.5)*GND_SIG
  CALL RUFFSOIL(FREQ(K)*1E9,GND_MV,GND_TEMP,KSIG,TETAD_SOIL,&
    EPS_UPPER,S0H,S0V)

  ! D.  COMPUTE SCATTERING COEFFICIENTS, REFLECTIVITIES AND 
  !     TRANSMISSIVITIES

  DEALLOCATE(RTEI,RDEI) ! BECAUSE THESE ARRAYS WILL CHANGE SIZES
  ALLOCATE(GBIH(RNUM),GBIV(RNUM),GS6(RNUM),GA2I(RNUM),&
    TSCAT(RNUM),RTEI(RNUM),RDEI(RNUM),RSIHLONG(RNUM+1),&
    RSIVLONG(RNUM+1))
     
  CALL SCCOEFF(RROI,RTI,RPCI,FREQ(K),RWIFR,RGAI,SCCHO,GBIH,GBIV,& 
    GS6,GA2I,RNUM,REPSI,REPSII,pixel,replicate,rank,meas)

  CALL PFADC(TETA,RDI,REPSI,GS6,RDEI,RTEI,TSCAT,RNUM)

  RSIHLONG(1)=S0H
  RSIHLONG(2:RNUM+1)=RSIH
  RSIVLONG(1)=S0V
  RSIVLONG(2:RNUM+1)=RSIV
  CALL POLMIX(TSCAT,RSIHLONG,RSIVLONG,RNUM)

  ! E.  COMPUTE BRIGHTNESS TEMPERATURES AND EMISSIVITIES OF SNOWPACK AND 
  !       GROUND

  ALLOCATE(RI(1:RNUM),TRI(1:RNUM),DH(1:RNUM),DV(1:RNUM))

  TSKYH=TB_UBC(1,K)
  TSKYV=TB_UBC(2,K)

  CALL RT(GA2I,GBIH,DEI,RI,TRI,RNUM)
  CALL LAYER(RI,RSIHLONG,TRI,RTI,GND_TEMP,TSKYH,DH,RNUM,pixel)

  TB_OUT(1,K)=(1-RSIHLONG(RNUM+1))*DH(RNUM)+RSIHLONG(RNUM+1)*TSKYH
  CALL RT(GA2I,GBIV,RDEI,RI,TRI,RNUM)
  CALL LAYER(RI,RSIVLONG,TRI,TI,GND_TEMP,TSKYV,DV,RNUM,pixel)
  TB_OUT(2,K)=(1-RSIVLONG(RNUM+1))*DV(RNUM)+RSIVLONG(RNUM+1)*TSKYV
! CALL EMISSIVITY(GA2I,GBIH,GBIV,RDEI,RSIHLONG,RSIVLONG,GND_TEMP,RTI,&
!   ESG_H,ESG_V,RNUM)
 
  ! F.  DEALLOCATE VARIABLE SPACE

  DEALLOCATE(ROI,GDIMM,PCI,EPSI,& 
    EPSII,GAI,NS,TEI,DEI,SIH,SIV,EPSR,GBIH,GBIV,GS6,GA2I,TSCAT,&
    RSIHLONG,RSIVLONG,RI,TRI,DH,DV)

  DEALLOCATE(RROI,REPSI,REPSII,RTEI,RSIH,RSIV,RDI,RDEI,RTI,RPCI,RWIFR,RGAI,&
    XROI,XEPSI,XEPSII,XTEI,XSIH,XSIV,XDI,XDEI,XTI,XPCI,XWIFR,XGAI)

END DO

END SUBROUTINE SS_MODEL

! -------------------------------------------------------------------------
!
SUBROUTINE ABSCOEFF(EPSI,EPSII,TI,FREQ,WIFR,GAI,NUM)
!
! -------------------------------------------------------------------------
!
!   CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   COMPUTES THE ABSORPTION COEFFICIENT FROM THE DIELECTRIC PROPERTIES   
!
!
!   [GAI] = ABSCOEFF(EPSI,EPSII,TI,FREQ,WI)
!       GAI:   ABSORPTION COEFFICIENT [M^-1]
!       EPSI:  REAL PART DIEL
!       EPSII: IMAGINARY PART DIEL
!       TI:    PHYSICAL TEMPERATURE
!       FREQ:  FREQUENCY [GHZ]
!       WI:    VOLUMETRIC LIQUID WATER CONTENT
!
!   VERSION HISTORY:
!      1.0    WI 15.7.95
!      1.1    WI 12.11.97 MORE PRECISE FORMULA FOR GAI USED
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!
!   USES: NONE

IMPLICIT NONE

INTEGER, INTENT(IN) :: NUM
REAL, INTENT(IN) :: EPSI(NUM),EPSII(NUM),TI(NUM),FREQ,WIFR(NUM)
REAL, INTENT(OUT) :: GAI(NUM)
REAL C,PI

C=2.99793
PI=3.14159
GAI=((2*PI*10*FREQ)*EPSII)/(C*(EPSI-(EPSII**2/4*EPSI))**0.5)

END SUBROUTINE ABSCOEFF

! -------------------------------------------------------------------------
!
SUBROUTINE AOT(EPS_T,EPS_B,TETAD_SOIL)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM PULLIAINEN
!
!     FUNCTION FOR CALCULATING THE ANGLE OF TRANSMISSION 
!     IN THE CASE OF A PLANE BOUNDARY BETWEEN TWO
!     MEDIAS WITH COMPLEX PERMITTIVITIES
!     K. TIGERSTEDT/97
!
!   INPUT PARAMETERS:
!     EPS_T = TOP LAYER EPSILON_R [COMPLEX],
!     EPS_B = BOTTOM LAYER EPSILON_R [COMPLEX]
!     THETAD = INCIDENCE ANGLE [DEG]
!   OUTPUT:
!     KSID   = TRANSMISSION ANGLE IN BOTTOM LAYER [DEG]
!
!   VERSION HISTORY:
!      1.0    KT ?.?.97
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB

IMPLICIT NONE

REAL, INTENT(IN) :: EPS_T,EPS_B
REAL, INTENT(INOUT) :: TETAD_SOIL
REAL THETA,KSI,PI
COMPLEX NUMERATOR,DENOMINATOR

PI=3.14159
THETA=TETAD_SOIL/180*PI
NUMERATOR=EPS_T**0.5*SIN(THETA)
DENOMINATOR=(EPS_B-EPS_T*(1-COS(THETA)**2))**0.5

KSI=ATAN(REAL(NUMERATOR)/REAL(DENOMINATOR))
TETAD_SOIL=KSI/PI*180

END SUBROUTINE AOT

! -------------------------------------------------------------------------
!
SUBROUTINE BORN(K,VFI,PCIMM,EPSI,EICE,GB6,GC6,GF6,GS6,NUM,pixel,replicate,&
   rank,meas)
!
! -------------------------------------------------------------------------
!
!   CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER, MODIFIED BY MIKE... SEE 
!     VERSION HISTORY
!
!   CALCULATES THE SCATTERING COEFFICIENT USING BORN APPROXIMATION       
!
!   [GB6,GC6,GF6,GS6] = BORNA(K,VFI,PCI,EPSI,EICE,EPSEFF,KP)
!       GB6: 6-FLUX BACK SCATTERING COEFFICIENT 
!       GC6: 6-FLUX CROSS SCATTERING COEFFICIENT 
!       GF6: 6-FLUX FORWARD SCATTERING COEFFICIENT 
!       GS6: 6-FLUX SCATTERING COEFFICIENT 
!       K:   WAVE NUMBER
!       VFI: VOLUME FRACTION OF ICE
!       PCI: CORRELATION LENGTH
!       EPSI: DIELECTRIC CONSTANT OF SNOW
!       EICE: DIELECTRIC CONSTANT OF SNOW
!       EPSEFF: EFFECTIVE PERMITTIVITY
!
!   VERSION HISTORY:
!      1.0     WI 27.05.98 
!      2.0     WI 31.03.05.  TRANSLATED TO FORTRAN BY MIKE.  BORNA AND 
!                             BORNSNK COMBINED INTO ONE.  CALL TO POLDER
!                             IS USED TO COMPUTE EFFECTIVE PERMITTIVITY
!      3.0     MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!
!   USES:
!     
!   COPYRIGHT (C) 1998 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER, INTENT(IN) :: NUM
REAL, INTENT(IN) :: K,VFI(NUM),PCIMM(NUM),EPSI(NUM),EICE
REAL, INTENT(OUT) :: GB6(NUM),GC6(NUM),GF6(NUM),GS6(NUM)
REAL,DIMENSION(:),ALLOCATABLE :: PCI,A,A3,EA,EA3,K1,K3,KP,MUC,AA,XX,BB,BT,BF,&
  BTOT, EPSEFF,NEGMUC
INTEGER STEPS,ROW,ARG_LENGTH(2),I

integer,intent(in) :: pixel,replicate,rank,meas

ALLOCATE( PCI(NUM),A(NUM),A3(NUM),EA(NUM),EA3(NUM),K1(NUM),K3(NUM),KP(NUM),&
       MUC(NUM),AA(NUM),XX(NUM),BB(NUM),BT(NUM),BF(NUM),BTOT(NUM),&
       EPSEFF(NUM),NEGMUC(NUM) )
  
! 0)  CONSTANTS AND CONVERSION
STEPS=11
PCI=PCIMM*0.001

! 1)  COMPUTE FIELD FACTOR AND DEPOLARIZATION RATIO, FROM BORNSNK
!     'A' AFTER NOTE 10, MATZLER 1997, COMMENT AND CALL FROM BORNSNK
CALL SNOWAO(VFI,A,NUM)

! 2)  COMPUTE EFFECTIVE PERMITTIVITY USING NEW FUNCTION POLDER.M

CALL POLDER(VFI,A,EICE,EPSI,EPSEFF,NUM)

! 3)  COMPUTE KP, THE SQUARED RATIO BETWEEN INTERNAL / EXTERNAL FIELDS
!     THIS CODE ORIGINALLY IN BORNSNK, STARTING AT LINE 41
A3=1-2*A
EA=EPSEFF*(1-A)+A
EA3=EPSEFF*(1-A3)+A3
K1=(EA/(EA+A*(EICE-1)))**2
K3=(EA3/(EA3+A3*(EICE-1)))**2
KP=(2*K1+K3)/3

! 4)  COMPUTE SCATTERING COEFFICIENTS, FROM BORNA, STARTING AT LINE 34

!NOTE: THE NEXT LINE CONSISTENT WITH MATZLER AND WIESMANN 99 PAPER EQUATION (7).
!  BASED ON CORRESPONDENCE WITH MATZLER, EPSEFF=N^2

MUC=((EPSEFF-1)/EPSEFF)**0.5      

AA=2*(PCI*K)**3*K*VFI*(1-VFI)*(EICE-1)**2*KP
XX=PCI*K*EPSI**0.5

! 5)  TRIPLE INTEGRATION 

!  A) BACKWARD SCATTERING
ARG_LENGTH=(/1, NUM/)
NEGMUC=-1.*MUC
!CALL INTEGRMUI(XX,MUC,1.,-1.,-1.*MUC,STEPS,NUM,BB,ARG_LENGTH)
CALL INTEGRMUI(XX,MUC,1.,-1.,NEGMUC,STEPS,NUM,BB,ARG_LENGTH)
!  B) TRANSVERSE SCATTERING 
ARG_LENGTH=(/NUM, NUM/)
!CALL INTEGRMUI(XX,MUC,1.,-1.*MUC,MUC,STEPS,NUM,BT,ARG_LENGTH)
CALL INTEGRMUI(XX,MUC,1.,NEGMUC,MUC,STEPS,NUM,BT,ARG_LENGTH)
!  C) FORWARD SCATTERING   
ARG_LENGTH=(/NUM, 1/)
CALL INTEGRMUI(XX,MUC,1.,MUC,1.,STEPS,NUM,BF,ARG_LENGTH)

! 6)  CALCULATION OF SCATTERING COEFFICIENTS

BTOT=BB+BT+BF
GB6=AA*BB
GC6=0.25*AA*BT
GF6=AA*BF
GS6=AA*BTOT

DEALLOCATE( PCI,A,A3,EA,EA3,K1,K3,KP,MUC,AA,XX,BB,BT,BF,BTOT, EPSEFF,NEGMUC)

END SUBROUTINE BORN


! -------------------------------------------------------------------------
!
SUBROUTINE EMISSIVITY(GA2I,GBIH,GBIV,DEI,SIHLONG,SIVLONG,GND_TEMP,TI,&
                        EH,EV,NUM)
!
! -------------------------------------------------------------------------
!
!  CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER, MODIFIED BY MIKE. 
!    SEE VERSION HISTORY
!
!   CALCULATES THE SCATTERING COEFFICIENT USING BORN APPROXIMATION       
!
!       GA2I: ABSORPTION COEFFICIENT
!       GBIH: 2-FLUX SCATTERING COEFFICIENT, H POLARIZATION
!       GBIV: 2-FLUX SCATTERING COEFFICIENT, V POLARIZATION
!       DEI:  EFFECTIVE PATH LENGTH [M]
!       SIHLONG: LAYER INTERFACE REFLECTIVITY, H POLARIZATION
!       SIVLONG: LAYER INTERFACE REFLECTIVITY, V POLARIZATION
!       EH:   EMISSIVITY, H POLARIZATION
!       EV:   EMISSIVITY, V POLARIZATION
!
!   VERSION HISTORY:
!      1.0     MD 1 APR 05 THIS CODE WAS PART OF LMAIN.  I TRANSLATED TO 
!                 FORTRAN FROM MATLAB AND MOVED IT TO A SEPARATE SUBROUTINE.
!                 COMPARE WIESMANN AND MATZLER, 99 EQN (8)
!
!   USES: RT, LAYER
!     
!   COPYRIGHT (C) 1998 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER, INTENT(IN) :: NUM
REAL,INTENT(IN) :: GA2I(NUM),GBIH(NUM),GBIV(NUM),DEI(NUM),&
                     SIHLONG(NUM+1),SIVLONG(NUM+1),GND_TEMP,&
                     TI(NUM)
REAL,INTENT(OUT) :: EH,EV
REAL RI(NUM),TRI(NUM),DH(NUM),DV(NUM),TBH0,TBH100,TBV0,TBV100,TSKY

! HORIZONTAL BRIGHTNESS TEMPERATURES UNDER DIFFERENT TSKY VALS

CALL RT(GA2I,GBIH,DEI,RI,TRI,NUM)

TSKY=0.
CALL LAYER(RI,SIHLONG,TRI,TI,GND_TEMP,TSKY,DH,NUM)
TBH0=(1-SIHLONG(NUM+1))*DH(NUM)+SIHLONG(NUM+1)*TSKY

TSKY=100.
CALL LAYER(RI,SIHLONG,TRI,TI,GND_TEMP,TSKY,DH,NUM)
TBH100=(1-SIHLONG(NUM+1))*DH(NUM)+SIHLONG(NUM+1)*TSKY

! VERTICAL BRIGHTNESS TEMPERATUERS UNDER DIFFERENT TSKY VALS

CALL RT(GA2I,GBIV,DEI,RI,TRI,NUM)

TSKY=0.
CALL LAYER(RI,SIVLONG,TRI,TI,GND_TEMP,TSKY,DV,NUM)
TBV0=(1-SIVLONG(NUM+1))*DV(NUM)+SIVLONG(NUM+1)*TSKY

TSKY=100.
CALL LAYER(RI,SIVLONG,TRI,TI,GND_TEMP,TSKY,DV,NUM)
TBV100=(1-SIVLONG(NUM+1))*DV(NUM)+SIVLONG(NUM+1)*TSKY

! COMPUTE EMISSIVITIES
EH=1-(TBH100-TBH0)/100
EV=1-(TBV100-TBV0)/100

END SUBROUTINE EMISSIVITY

! -------------------------------------------------------------------------
!
SUBROUTINE EPSICE(TI,FREQ,EICE,NUM)
!
! -------------------------------------------------------------------------
!
!   CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   CALCULATES THE DIELECTRIC PERMITTIVITY OF ICE
!   AFTER HUFFORD, MITZIMA AND MATZLER
!
!   EICE = EPSICE(TI,FREQ)
!      EICE:  DIELECTRIC PERMITTIVITY OF ICE
!      TI:    TEMPERATURE IN K 
!      FREQ:  FREQUENCY IN GHZ
!
!   VERSION HISTORY:
!      1.0    WI 15.7.95
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!   
!   USES: NONE
!
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER, INTENT(IN) :: NUM
REAL, INTENT(IN) :: TI(NUM),FREQ
REAL, INTENT(OUT) ::  EICE(NUM)
REAL PP(NUM),B,B1,B2,DB(NUM),BETA(NUM),ALPHA(NUM)

PP=(300/TI)-1
B=335.25
B1=0.0207
B2=1.16E-11
DB=EXP(-10.02+0.0364*(TI-273))
BETA=((B1*EXP(B/TI))/(TI*(EXP(B/TI)-1)**2))+B2*FREQ**2+DB
ALPHA=(0.00504+0.0062*PP)*EXP(-22.1*PP)
EICE=ALPHA/FREQ+BETA*FREQ

END SUBROUTINE EPSICE

! -------------------------------------------------------------------------
!
SUBROUTINE EPSR(ROI,NUM,EPSI)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   CALCULATES THE DIELECTRIC PERMITTIVITY FOR DRY SNOW FROM 
!   DENSITY .
!
!   EPSI = EPSR(ROI)
!       EPSI:  REAL PART OF DIELECTRIC PERMITTIVITY
!       ROI:   DENSITY G/CM^3
!
!   VERSION HISTORY:
!      1.0    WI 15.7.95
!      1.1    WI 23.9.97 ADDED LOOYENGA FOR SNOW DENSER THAN 0.4 G/CM^3
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!
!   USES:
!       EPSICE
!
!
!
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER, INTENT(IN) :: NUM
REAL, INTENT(IN) :: ROI(NUM)
REAL, INTENT(OUT) :: EPSI(NUM)
INTEGER :: I
REAL :: EHB, ESB
REAL,DIMENSION(:),ALLOCATABLE :: VFI

ALLOCATE(VFI(NUM))

VFI=ROI/0.917
EHB=0.99913
ESB=1.4759

DO I=1,NUM
  IF (ROI(I)<=0.4) THEN
    EPSI(I)=1+1.5995*ROI(I)+1.861*ROI(I)**3
  ELSE
    EPSI(I)=((1-VFI(I))*EHB+VFI(I)*ESB)**3
  END IF
END DO   

DEALLOCATE(VFI)

END SUBROUTINE EPSR

! -------------------------------------------------------------------------
!
SUBROUTINE EPSS(MV,T,F,EPSS_VAR)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM PULLIAINEN
!
!     FUNCTION FOR CALCULATING EPSILON FOR SOIL USING FREQUENCY,
!     TEMPERATURE AND VOLUMETRIC SOIL MOSITURE. USES EPSW.M FOR
!     DIELECTRICITY OF WATER. 
!
!     BY J. PULLIAINEN (MOD. BY K. TIGERSTEDT) 
! 
!     MV [0..1]
!     T [C]
!     F [HZ]
!
!     NOTE: SOME FINNISH COMMENTS WERE NOT COPIED IN ENTIRETY -MD
!   VERSION HISTORY:
!      1.0    JP ?.?.?
!      1.1    KT ?.?.?
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB

IMPLICIT NONE

REAL,INTENT(IN) :: MV,T,F
COMPLEX,INTENT(OUT) :: EPSS_VAR
REAL S,CLAY,RHOS,BETA,ALFA,REI,A,B,C,M,IEI_S,A_P,B_P,C_P,IEI_P,&
      DELTA_IEI,INVT,ALF,B_1,B_2,BB,BET_M,BET_DELTA,BET,IEI,SS,EW_R,EW_I
COMPLEX EW,EPSALF

S=0.4                          ! SAND FRACTION
CLAY=0.3                       ! CLAY FRACTION
RHOS=1.1                       ! DRY-MATTER DENSITY

BETA=1.09-0.11*S+0.18*CLAY
ALFA=0.65

IF(T>0)THEN
  CALL EPSW(F,T,EW_R,EW_I)
ELSE
  REI=3.1884+9.1E-4*T          ! MATZLER AND WEGMULLER 1987
  A=0.0026                     ! IMPURE ICE -5 ASTETTA (MATZLER)
  B=0.00023
  C=0.87
  M=F/1E9
  IEI_S=A/M+B*M**C
  A_P=6E-4                     ! PURE ICE -5 ASTETTA (MATZLER)
  B_P=6.5E-5
  C_P=0.7
  IEI_P=A_P/M+B_P*M**C_P
  DELTA_IEI=IEI_S-IEI_P

  ! HUFFORD 1991        
  INVT=300/(T+273)-1
  ALF=(0.00504+0.0062*INVT)*EXP(-22.1*INVT)

  ! (MISHIMA,MATZLER)
  B_1=0.0207
  B_2=1.16E-11
  BB=335
  BET_M=B_1/(T+273)*EXP(BB/(T+273))/(EXP(BB/(T+273))-1)**2+B_2*M**2
  BET_DELTA=EXP(-10.02+0.0364*T)
  BET=BET_M+BET_DELTA
  IEI=ALF/M+BET*M
  SS=10
  IEI=IEI+DELTA_IEI*SS/13
  EW_R=REI
  EW_I=IEI
END IF

EW=CMPLX(EW_R,(-1*EW_I))

EPSALF=1+0.65*RHOS+MV**BETA*(EW**ALFA-1)

EPSS_VAR=EPSALF**(1/ALFA) 

END SUBROUTINE EPSS


! -------------------------------------------------------------------------
!
SUBROUTINE EPSW(F,T,EW_R,EW_I)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM PULLIAINEN
!
!     LASKEE PUHTAAN VEDEN EPSILONIN TAAJUUDEN F, 
!     L MP TILAN T FUNKTIONA
!     5.12.89 JPK; 17.2.1993 JP
!    VERSION HISTORY:
!      1.0    JPK 5.12.89 
!      1.1    JP 17.02.93
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!
!     NOTE: ADDITIONAL COMMENTS IN FINNISH NOT COPIED HERE -MD

IMPLICIT NONE

REAL,INTENT(IN) :: F,T
REAL,INTENT(OUT) :: EW_R,EW_I
REAL S,E0,EW_INF,EW0,N,D,ALFA,SIGMA_25,SIGMA_VAR,TW,PI,SIGMA

PI=3.14159
S=0
E0=8.854E12
EW_INF=4.9
EW0=87.74-0.40008*T+9.398E-4*T**2+1.410E-6*T**3
N=0

D=25-T
ALFA=2.033E-2+1.266E-4*D+2.462E-6*D**2-S*(1.849E-5-2.551E-7*D+&
       2.551E-8*D**2)

SIGMA_25=S*(0.182521-1.46192E-3*S+2.09324E-5*S**2-1.28205E-7*S**3)
SIGMA=SIGMA_25*EXP(-D*ALFA)

TW=1/(2*PI)*(1.1109E-10-3.824E-12*T+6.938E-14*T*T-5.096E-16*T*T*T)

EW_R=EW_INF+(EW0-EW_INF)/(1+(2*PI*F*TW)**2)
EW_I=(EW0-EW_INF)*2*PI*F*TW/(1+(2*PI*F*TW)**2)+SIGMA/(2*PI*E0*F)

END SUBROUTINE EPSW

! -------------------------------------------------------------------------
!
SUBROUTINE FRESNELC(TEI,EPSI,SIH,SIV,NUM)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   FRESNEL REFLECTION COEFFICIENTS (ASSUMING EPS'' = 0)
!     (LAYER N+1 IS THE AIR ABOVE THE SNOWPACK)
!
!   [SIH,SIV] = FRESNEL(TEI,ROI)
!       SIH:  INTERFACE REFLECTIVITY AT H POL
!       SIV:  INTERFACE REFLECTIVITY AT V POL
!       TEI:  LOCAL INCIDENCE ANGLE
!       EPSI: REAL PART OF DIELECTRIC PERMITTIVITY
!
!   VERSION HISTORY:
!      1.0    WI 15.7.97
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!   
!   USES:
!       EPSR
!
!
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER, INTENT(IN) :: NUM
REAL, INTENT(IN) :: TEI(NUM+1),EPSI(NUM)
!ctvs REAL, INTENT(OUT) :: SIH(NUM-1),SIV(NUM-1)
 REAL, INTENT(OUT) :: SIH(NUM),SIV(NUM)
INTEGER :: I
REAL EPSO,EPSU,TEIN
REAL,DIMENSION(:),ALLOCATABLE:: EPSI_LOCAL

ALLOCATE(EPSI_LOCAL(NUM+1))

EPSI_LOCAL(1:NUM)=EPSI
EPSI_LOCAL(NUM+1)=1

DO I=1,NUM
  EPSO=EPSI_LOCAL(I+1)
  EPSU=EPSI_LOCAL(I)
  TEIN=TEI(I+1)
  SIH(I)=((EPSO**0.5*COS(TEIN)-(EPSU-EPSO*SIN(TEIN)**2)**0.5)/&
          (EPSO**0.5*COS(TEIN)+(EPSU-EPSO*SIN(TEIN)**2)**0.5))**2
  SIV(I)=((EPSU*COS(TEIN)-EPSO**0.5*(EPSU-EPSO*SIN(TEIN)**2)**0.5)/&
          (EPSU*COS(TEIN)+EPSO**0.5*(EPSU-EPSO*SIN(TEIN)**2)**0.5))**2
END DO        

DEALLOCATE(EPSI_LOCAL)

END SUBROUTINE FRESNELC

! -------------------------------------------------------------------------
!
SUBROUTINE FRESNELRC(TEI,EPSI,FH,FV,NUM)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!
!   FRESNEL REFLECTION COEFFICIENTS (ASSUMING EPS'' = 0)
!     (LAYER N+1 IS THE AIR ABOVE THE SNOWPACK)
!
!   [FH,FV] = FRESNELRC(TEI,EPSR)
!       FH:   FRESNEL REFLECTION COEFFICIENT AT H POL
!       FV:   FRESNEL REFLECTION COEFFICIENT AT V POL
!       TEI:  LOCAL INCIDENCE ANGLE
!       EPSR: (REAL PART) DIELECTRIC PERMITTIVITY
!
!   VERSION HISTORY:
!      1.0    WI 15.7.95
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!   
!   USES:
!
!
!
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER,INTENT(IN) :: NUM
REAL,INTENT(IN) :: TEI(NUM+1),EPSI(NUM+1)
REAL,INTENT(OUT) :: FH(NUM),FV(NUM)
INTEGER N
REAL EPSN,TEIN,SINQ,QEPS,WURZ,WSUB,ND

DO N=1,NUM
  EPSN=EPSI(N)/EPSI(N+1)
  TEIN=TEI(N+1)
  SINQ=SIN(TEIN)**2
  QEPS=SINQ/EPSN
  WURZ=(1-QEPS)**0.5
  WSUB=EPSN-SINQ
  ND=EPSN**0.5
  
  FH(N)=((ND*WURZ-COS(TEIN))/(ND*WURZ+COS(TEIN)))
  FV(N)=((WURZ-ND*COS(TEIN))/(WURZ+ND*COS(TEIN)))
END DO             

END SUBROUTINE FRESNELRC

! -------------------------------------------------------------------------
!
SUBROUTINE GAMMAH(EPSS,THETA,GAMMAH_VAR)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM PULLIAINEN
!
!    VERSION HISTORY:
!      1.0    JPK 5.12.89 
!      1.1    JP 17.03.93
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!
!     NOTE: ADDITIONAL FINNISH COMMENTS NOT COPIED HERE

IMPLICIT NONE

REAL,INTENT(IN) :: THETA
COMPLEX,INTENT(IN) :: EPSS      
REAL,INTENT(OUT) :: GAMMAH_VAR
REAL THETA_RAD,PI,COSTHETA
COMPLEX NELIO

PI=3.14159
THETA_RAD=THETA/180*PI
COSTHETA=COS(THETA_RAD)
NELIO=SQRT(EPSS-SIN(THETA_RAD)**2)
GAMMAH_VAR=(ABS((COSTHETA-NELIO)/(COSTHETA+NELIO)))**2

END SUBROUTINE GAMMAH

! -------------------------------------------------------------------------
!
SUBROUTINE GAMMAV(EPSS,THETA,GAMMAV_VAR)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM PULLIAINEN
!
!    VERSION HISTORY:
!      1.0    JPK 5.12.89 
!      1.1    JP 17.03.93
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!
!     NOTE: ADDITIONAL FINNISH COMMENTS NOT COPIED HERE

IMPLICIT NONE

REAL,INTENT(IN) :: THETA
COMPLEX,INTENT(IN) :: EPSS      
REAL,INTENT(OUT) :: GAMMAV_VAR
REAL THETA_RAD,PI,COSTHETA
COMPLEX NELIO

PI=3.14159
THETA_RAD=THETA/180*PI
COSTHETA=COS(THETA_RAD)
NELIO=SQRT(EPSS-SIN(THETA_RAD)**2)
GAMMAV_VAR=(ABS((EPSS*COSTHETA-NELIO)/(EPSS*COSTHETA+NELIO)))**2

END SUBROUTINE GAMMAV

! -------------------------------------------------------------------------
!
SUBROUTINE INTEGRFI(XX,MUI,MUO,STEPS,INTEGR,NUM)
!
! -------------------------------------------------------------------------
!
!   CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   CALCULATES INTEGRATION  OVER INCIDENT DIRECTIONS, MUI
!       FROM MINI TO MAXI IN STEPS INTERVALS OF THE FI
!       PHASE FUNCTION.
!   
!
!   INTEGR = INTEGRFI(XX,MUI,MUO,STEPS)
!
!   VERSION HISTORY:
!      1.0     WI 27.05.98 
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!      2.1    MD 21 NOV 05 CHANGED ALL LOCALS TO ALLOCATABLE
!
!   USES: - NONE
!     
!   COPYRIGHT (C) 1998 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER,INTENT(IN) :: NUM,STEPS
REAL, INTENT(IN) :: XX(NUM),MUI(NUM),MUO(NUM)
REAL, INTENT(OUT) :: INTEGR(NUM)
REAL PI,DELTA,FI,F0,COFI,SI2FI

INTEGER :: IFI

REAL,DIMENSION(:),ALLOCATABLE ::  X2,SII,SIO,COSTE,SI2CHI,FUNC

ALLOCATE( X2(NUM),SII(NUM),SIO(NUM),COSTE(NUM),SI2CHI(NUM),FUNC(NUM) )

PI=3.14159
DELTA=PI/STEPS
F0=0.5*DELTA
INTEGR=0
X2=2*XX**2

DO IFI=1,STEPS
  FI=F0+(IFI-1)*DELTA
  SII=(1-MUI**2)**0.5
  SIO=(1-MUO**2)**0.5
  COFI=COS(FI)
  SI2FI=1-COFI**2
  COSTE=MUI*MUO+SII*SIO*COFI
  SI2CHI=0.5*(1+COSTE**2)
  FUNC=SI2CHI/(1+(1-COSTE)*X2)**2
  INTEGR=INTEGR+FUNC*DELTA
END DO

INTEGR=INTEGR/PI

DEALLOCATE( X2,SII,SIO,COSTE,SI2CHI,FUNC )

END SUBROUTINE INTEGRFI

! -------------------------------------------------------------------------
!
SUBROUTINE INTEGRMUI(XX,MINI,MAXI,MINO,MAXO,STEPS,NUM,INTEGR,ARG_LENGTH)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   CALCULATES INTEGRATION  OVER INCIDENT DIRECTIONS, MUI
!       FROM MINI TO MAXI IN STEPS INTERVALS OF THE FI AND
!       MUO INTEGRATED PHASE FUNCTION.
!   
!
!   INTEGR = INTEGRMUI(XX,MINI,MAXI,MINO,MAXO,STEPS)
!
!   VERSION HISTORY:
!      1.0     WI 27.05.98 
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!
!   USES: INTEGRMUO
!     
!   COPYRIGHT (C) 1998 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE
 
INTEGER,INTENT(IN) :: NUM,STEPS,ARG_LENGTH(2)
REAL, INTENT(IN) :: XX(NUM),MINI(NUM),MAXI,MINO(ARG_LENGTH(1)),&
                      MAXO(ARG_LENGTH(2))
REAL, INTENT(OUT) :: INTEGR(NUM)
REAL :: STEPSR
REAL,DIMENSION(:),ALLOCATABLE :: DMU,DELTA,F0,MUI,FUNC
INTEGER IMU

ALLOCATE( DMU(NUM),DELTA(NUM),F0(NUM),MUI(NUM),FUNC(NUM) )

STEPSR=REAL(STEPS)
DMU=MAXI-MINI
DELTA=DMU/STEPSR
F0=0.5*DELTA
INTEGR=0.

DO IMU=1,STEPS
  MUI=MINI+F0+(IMU-1)*DELTA
  CALL INTEGRMUO(XX,MUI,MINO,MAXO,STEPS,FUNC,NUM,ARG_LENGTH)
  INTEGR=INTEGR+FUNC*DELTA
END DO

INTEGR=INTEGR/DMU

DEALLOCATE( DMU,DELTA,F0,MUI,FUNC )

END SUBROUTINE INTEGRMUI

! -------------------------------------------------------------------------
!
SUBROUTINE INTEGRMUO(XX,MUI,MINO,MAXO,STEPS,INTEGR,NUM,ARG_LENGTH)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   CALCULATES INTEGRATION  OVER INCIDENT DIRECTIONS, MUI
!       FROM MINO TO MAXO IN STEPS INTERVALS OF THE FI
!       PHASE FUNCTION.
! 
!
!   INTEGR = INTEGRMUO(XX,MUI,MINO,MAXO,STEPS)
!
!   VERSION HISTORY:
!      1.0     WI 27.05.98 
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!      2.1    MD 21 NOV 05 MADE ALL LOCALS ALLOCATABLE
!
!   USES: INTEGRFI
!     
!   COPYRIGHT (C) 1998 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND      

IMPLICIT NONE

INTEGER,INTENT(IN) :: NUM,STEPS,ARG_LENGTH(2)
REAL,INTENT(IN) :: MUI(NUM),MINO(ARG_LENGTH(1)),MAXO(ARG_LENGTH(2)),XX(NUM)
REAL,INTENT(OUT) :: INTEGR(NUM)
REAL :: MINOS,MAXOS
REAL,DIMENSION(:), ALLOCATABLE ::  DMU,DELTA,F0,MUO,FUNC
INTEGER IMU,I

ALLOCATE( DMU(NUM),DELTA(NUM),F0(NUM),MUO(NUM),FUNC(NUM) )

IF (ARG_LENGTH(1)==1) THEN
! N.B. IF MINO IS SCALAR, MAXO IS VECTOR
  MINOS=MINO(1)
  DMU=MAXO-MINOS
ELSEIF (ARG_LENGTH(2)==1) THEN
! N.B. IF MAXO IS SCALAR, MINO IS VECTOR
  MAXOS=MAXO(1)
  DMU=MAXOS-MINO
ELSE
  DMU=MAXO-MINO
END IF     

DELTA=DMU/STEPS
F0=0.5*DELTA
INTEGR=0.

DO IMU=1,STEPS
  IF (ARG_LENGTH(1)==1) THEN
    MUO=MINOS+F0+(IMU-1)*DELTA
  ELSE
    MUO=MINO+F0+(IMU-1)*DELTA
  END IF
  CALL INTEGRFI(XX,MUI,MUO,STEPS,FUNC,NUM)
  INTEGR=INTEGR+FUNC*DELTA
END DO

INTEGR=INTEGR/2

DEALLOCATE( DMU,DELTA,F0,MUO,FUNC )

END SUBROUTINE INTEGRMUO

! -------------------------------------------------------------------------
!
SUBROUTINE LAYER(RI,SI,TRI,TI,TGND,TSKY,D,NUM,pixel)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   CALCULATES THE UPWELLING BRIGHTNESS TEMPERATURES D (SEE NOTE 6)  
!
!   D = LAYER(RI,SI,TI,TI,TGND,TSKY)
!       D:    UPWELLING BRIGHTNESS TEMPERATURE
!       RI:   LAYER REFLECTIVITY
!       SI:   INTERFACE REFLECTIVITY
!       TI:   LAYER TRANSMISSIVITY
!       TI:   PHYSICAL TEMPERATURE [K]
!       TGND: BRIGHTNESS TEMPERATURE OF THE SOIL BELOW THE SNOWPACK
!       TSKY: BRIGHTNESS TEMPERATURE OF THE SKY
!
!   VERSION HISTORY:
!      1.0    WI 15.7.95
!      1.1    WI 26.9.97  HANDLES ALSO THE SPECIAL CASE OF A SINGLE LAYER NOW
!      1.2    WI 02.03.99 FIXED ERROR IN 1 LAYER HANDLING 
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!  
!   USES: - INVERT_MATRIX SUBROUTINE BY MIKE,ZEROSI INTERNAL FUNCTION AND
!              EYEI INTERNAL FUNCTION BY MIKE
!
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER,INTENT(IN) :: NUM
REAL,INTENT(IN) :: RI(NUM),SI(NUM+1),TRI(NUM),TI(NUM),TGND,TSKY
REAL,INTENT(OUT) :: D(NUM)
REAL K1
REAL,DIMENSION(:),ALLOCATABLE ::  EI,F,E
REAL,DIMENSION(:,:),ALLOCATABLE :: M1,H,EYE,M2,M3,M4,M5,INVIM1,INVIM5,EYEM1,&
  EYEM5 
INTEGER ROW,COL,I

integer,intent(in) :: pixel

ALLOCATE( EI(NUM),M1(NUM,NUM),H(NUM-1,NUM-1),EYE(NUM,NUM),M2(NUM,NUM),&
       M3(NUM,NUM),M4(NUM,NUM),E(NUM),F(NUM),M5(NUM,NUM),&
       INVIM1(NUM,NUM),INVIM5(NUM,NUM),EYEM1(NUM,NUM),EYEM5(NUM,NUM) )

EI=1-RI-TRI

IF (NUM==1) THEN
  !FOR ONE LAYER CASE, THIS IS THE ENTIRE CALCULATION; THE MATRIX INVERSION
  ! ROUTINE IS NOT NEEDED.  THEREFORE, A RETURN STATEMENT IS USED HERE...
  K1=(1-RI(1)*SI(1))*(1-RI(1)*SI(2))-TRI(1)*SI(1)*TRI(1)*SI(2)
  D=TRI(1)*SI(1)*((1-SI(1))*RI(1)*TGND+(1-SI(2))*TSKY*TRI(1)+&
      EI(1)*TI(1))/K1+(1-RI(1)*SI(1))*((1-SI(1))*TGND*TRI(1)+&
      (1-SI(2))*TSKY*RI(1)+EI(1)*TI(1))/K1
  RETURN
ELSE
! INITIALIZE ARRAYS
  M1=ZEROSI(NUM,NUM)
  H=ZEROSI(NUM-1,NUM-1)
  M2=ZEROSI(NUM,NUM)
  M3=ZEROSI(NUM,NUM)
  M4=ZEROSI(NUM,NUM)

  DO I=1,NUM
    M1(I,I)=RI(I)*SI(I)
  END DO
  DO I=1,NUM-1
    H(I,I)=TRI(I)*(1-SI(I+1))
  END DO
  M1(1:NUM-1:1,2:NUM:1)=M1(1:NUM-1:1,2:NUM:1)+H

  EYE=EYEI(NUM)
  H=ZEROSI(NUM-1,NUM-1)
  DO I=1,NUM
    M2(I,I)=TRI(I)*SI(I+1)
  END DO
  DO I=2,NUM
    H(I-1,I-1)=RI(I)*(1-SI(I))
  END DO
  M2(2:NUM:1,1:NUM-1:1)=M2(2:NUM:1,1:NUM-1:1)+H

  H=ZEROSI(NUM-1,NUM-1)
  DO I=1,NUM
    M3(I,I)=TRI(I)*SI(I)
  END DO
  DO I=1,NUM-1
    H(I,I)=RI(I)*(1-SI(I+1))
  END DO
  M3(1:NUM-1:1,2:NUM:1)=M3(1:NUM-1:1,2:NUM:1)+H

  H=ZEROSI(NUM-1,NUM-1)
  DO I=1,NUM
    M4(I,I)=RI(I)*SI(I+1)
  END DO
  DO I=2,NUM
    H(I-1,I-1)=TRI(I)*(1-SI(I))
  END DO
  M4(2:NUM:1,1:NUM-1:1)=M4(2:NUM:1,1:NUM-1:1)+H

  E=EI*TI
  E(1)=E(1)+RI(1)*(1-SI(1))*TGND
  E(NUM)=E(NUM)+TRI(NUM)*(1-SI(NUM+1))*TSKY

  F=EI*TI
  F(1)=F(1)+TRI(1)*(1-SI(1))*TGND
  F(NUM)=F(NUM)+RI(NUM)*(1-SI(NUM+1))*TSKY

END IF
   
EYEM1=EYE-M1
CALL INVERT_MATRIX(EYEM1,INVIM1,NUM)

M5=MATMUL(M3,MATMUL(INVIM1,M2))+M4
EYEM5=EYE-M5
CALL INVERT_MATRIX(EYEM5,INVIM5,NUM)

D=MATMUL(INVIM5,(MATMUL(MATMUL(M3,INVIM1),E)+F))

DEALLOCATE( EI,M1,H,EYE,M2,M3,M4,E,F,M5,INVIM1,INVIM5,EYEM1,EYEM5 )

CONTAINS

FUNCTION ZEROSI(ROW,COL)

! BY MIKE, 1 APRIL 2005
!
! THIS FUNCTION IS INTERNAL TO LAYER SUBROUTINE
! MANY VALUES IN THE ABOVE ARRAYS ARE NEVER A VALUE OTHER THAN ZERO, SO INSTEAD
! OF USING WHATEVER VALUES ARE RANDOMLY ALLOCATED TO THE ARRAYS WHEN THEY ARE 
! ALLOCATED, I SPECIFICALLY SET EACH POSITION TO ZERO.  THIS FUNCTION SHOULD 
! BE ENTIRELY unneceSSARY, SINCE THE DEFAULT VALUES IN THE ARRAY SHOULD BE ZERO. NONETHELESS, I DECIDED TO EXPLICITLY ZERO OUT 2-D ARRAYS BEFORE USING THEM

IMPLICIT NONE

INTEGER,INTENT(IN):: ROW,COL
REAL :: ZEROSI(ROW,COL)
INTEGER I,J

DO I=1,ROW
  DO J=1,COL
    ZEROSI(I,J)=0.
  END DO
END DO

END FUNCTION ZEROSI

FUNCTION EYEI(NUM)


! BY MIKE, 1 APRIL 2005
!
! THIS FUNCTION IS INTERNAL TO LAYER SUBROUTINE
! THIS FUNCTION DEFINES AN IDENTITY MATRIX OF SIZE NUM

INTEGER, INTENT(IN) :: NUM
REAL :: EYEI(NUM,NUM)
INTEGER I,J
DO I=1,NUM
  DO J=1,NUM
    IF (I==J) THEN
      EYEI(I,J)=1
    ELSE
      EYEI(I,J)=0
    END IF
  END DO
END DO
END FUNCTION EYEI
  

END SUBROUTINE LAYER

! -------------------------------------------------------------------------
!
SUBROUTINE MIXMOD(FREQ,TI,WIFR,EPSI,EPSII,NUM)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   CALCULATES THE PERMITTIVITY FOR WETNESS > 0
!      PHYSICAL MIXING MODEL WEISE 97 AFTER MATZLER 1987 (CORRECTED)
!      WATER TEMPERATURE IS ASSUMED CONSTANT AT 273.15 K
!
!   [EPSI,EPSII] = MIXMOD(F,TI,WI,EPSI,EPSII)
!       EPSI:  REAL PART OF THE PERMITTIVITY
!       EPSII: IMAGINARY PART OF THE PERMITTIVITY
!       F:     FREQUENCY [GHZ]
!       TI:    PHYSICAL SNOW TEMPERATURE
!       WI:    WETNESS [!]
!       EPSI:  REAL PART OF DRY SNOW PERM.
!       EPSII: IMAGINARY PART OF DRY SNOW PERM.
!
!   VERSION HISTORY:
!      1.0    WI 15.7.95
!      2.0    MD  1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!      2.1    MD 21 NOV 05 MADE ALL LOCALS DYNAMICALLY ALLOCATABLE
!   
!   USES: - NONE
!
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER, INTENT(IN) :: NUM
REAL, INTENT(IN) :: FREQ, TI(NUM), WIFR(NUM)
REAL, INTENT(INOUT) :: EPSI(NUM), EPSII(NUM)

REAL :: AA,AB,AC,EUW,ESW,FRW
REAL,DIMENSION(:),ALLOCATABLE ::  ESA,ESB,ESC,EUA,EUB,EUC,FA,FB,FC,EEA,EEB,&
     EEC,FWA,FWB,FWC,DEPSIA,DEPSIB,DEPSIC,DEPSI,DEPSIIA,DEPSIIB,DEPSIIC,DEPSII

ALLOCATE( ESA(NUM),ESB(NUM),ESC(NUM),EUA(NUM),&
     EUB(NUM),EUC(NUM),FA(NUM),FB(NUM),FC(NUM),EEA(NUM),EEB(NUM),&
     EEC(NUM),FWA(NUM),FWB(NUM),FWC(NUM),DEPSIA(NUM),DEPSIB(NUM),&
     DEPSIC(NUM),DEPSI(NUM),DEPSIIA(NUM),DEPSIIB(NUM),DEPSIIC(NUM),&
     DEPSII(NUM) )

AA=0.005
AB=0.4975
AC=0.4975
EUW=4.9
ESW=88.045
FRW=0.11109                            ! INVERSE RELAXATION FREQUENCY OF WATER
ESA=(ESW-EPSI)/(3*(1+AA*(ESW/EPSI-1)))
ESB=(ESW-EPSI)/(3*(1+AB*(ESW/EPSI-1)))
ESC=(ESW-EPSI)/(3*(1+AC*(ESW/EPSI-1)))
EUA=(EUW-EPSI)/(3*(1+AA*(EUW/EPSI-1)))
EUB=(EUW-EPSI)/(3*(1+AB*(EUW/EPSI-1)))
EUC=(EUW-EPSI)/(3*(1+AC*(EUW/EPSI-1)))

FA=1+AA*(ESW-EUW)/(EPSI+AA*(EUW-EPSI))
FB=1+AB*(ESW-EUW)/(EPSI+AB*(EUW-EPSI))
FC=1+AC*(ESW-EUW)/(EPSI+AC*(EUW-EPSI))

EEA=ESA-EUA
EEB=ESB-EUB
EEC=ESC-EUC

FWA=FRW/FA
FWB=FRW/FB
FWC=FRW/FC

DEPSIA=EUA+EEA/(1+(FWA*FREQ)**2)
DEPSIB=EUB+EEB/(1+(FWB*FREQ)**2)
DEPSIC=EUC+EEC/(1+(FWC*FREQ)**2)
DEPSI=WIFR*(DEPSIA+DEPSIB+DEPSIC)

DEPSIIA=FWA*FREQ*EEA/(1+(FWA*FREQ)**2)
DEPSIIB=FWB*FREQ*EEB/(1+(FWB*FREQ)**2)
DEPSIIC=FWC*FREQ*EEC/(1+(FWC*FREQ)**2)
DEPSII=WIFR*(DEPSIIA+DEPSIIB+DEPSIIC)

EPSI=EPSI+DEPSI
EPSII=EPSII+DEPSII

DEALLOCATE( ESA,ESB,ESC,EUA, EUB,EUC,FA,FB,FC,EEA,EEB,EEC,FWA,FWB,&
     FWC,DEPSIA,DEPSIB, DEPSIC,DEPSI,DEPSIIA,DEPSIIB,DEPSIIC, DEPSII )

END SUBROUTINE MIXMOD

! -------------------------------------------------------------------------
!
SUBROUTINE PFADC(TETA,DI,EPSI,GS6,DEI,TEI,TSCAT,NUM)
!
! -------------------------------------------------------------------------
!
!      CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   CALCULATES THE EFFECTIVE PATH LENGTH IN A LAYER
!
!   [DEI,TEI,TSCAT] = PFADC(TETA,DI,EPSI,GS6)
!       DEI:  EFFECTIVE PATH LENGTH [M]
!       TEI:  LOCAL INCIDENCE ANGLE
!       TSCAT: SCATTERING 
!       TETA: INCIDENCE ANGLE AT SNOW AIR INTERFACE
!       DI:   THICKNESS [M]
!       EPSI: DIELECTRIC PERMITTIVITY
!       GS6:  6-FLUX SCATTERING COEFFICIENT
!
!   VERSION HISTORY:
!      1.0    WI 15.10.97
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!      2.1    MD 21 NOV 05 CHANGED ALL LOCALS TO ALLOCATBLE
!
!   
!   USES: - NONE
!
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER, INTENT(IN) :: NUM
REAL, INTENT(IN) :: TETA,DI(NUM),EPSI(NUM),GS6(NUM)
REAL, INTENT(OUT) :: DEI(NUM),TEI(NUM),TSCAT(NUM)
REAL,DIMENSION(:),ALLOCATABLE ::  NS,COSTETASN,COSC,COSTETASC,TAUSCAT,COSTETA
INTEGER I

ALLOCATE( NS(NUM),COSTETASN(NUM),COSC(NUM),COSTETASC(NUM),TAUSCAT(NUM+1),&
       COSTETA(NUM) )

NS=EPSI**0.5
COSTETASN=(1-(SIN(TETA)/NS)**2)**0.5      
COSC=(1-(1/NS)**2)**0.5
COSTETASC=0.5*(1+COSC)
DEI=DI/COSTETASN
    
TAUSCAT(NUM+1)=0

DO I=NUM,1,-1
  TAUSCAT(I)=TAUSCAT(I+1)+DEI(I)*GS6(I)/2
  TSCAT(I)=EXP(-1*TAUSCAT(I))
  COSTETA(I)=TSCAT(I)*COSTETASN(I)+(1-TSCAT(I))*COSTETASC(I)
END DO      

TEI=ACOS(COSTETA)

DEALLOCATE (  NS,COSTETASN,COSC,COSTETASC,TAUSCAT,COSTETA )

END SUBROUTINE PFADC

! -------------------------------------------------------------------------
!
SUBROUTINE PFADI(TEI,DI,DEI,NUM)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   CALCULATES THE EFFECTIVE PATH LENGTH IN A LAYER
!
!   DEI = PFADI(TEI,DI)
!       DEI:  EFFECTIVE PATH LENGTH [M]
!       TEI:  LOCAL INCIDENCE ANGLE
!       DI:   THICKNESS [M]
!
!   VERSION HISTORY:
!      1.0    WI 15.7.95
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!   
!   USES - NONE
!
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER, INTENT(IN) :: NUM
REAL, INTENT(IN) :: TEI(NUM+1),DI(NUM)
REAL, INTENT(OUT) :: DEI(NUM)

DEI=DI/COS(TEI(1:NUM))

END SUBROUTINE PFADI

! -------------------------------------------------------------------------
!
SUBROUTINE POLDER(VFI,A,EICE,EPSI,EPSEFF,NUM)
!
! -------------------------------------------------------------------------
!  WRITTEN BY MIKE, 1 APRIL 05
!  SOLVES (1) THROUGH (5) IN MATZLER, 1996 IN IEEE
!  IMPLEMENTS NEWTON-RAPHSON ITERATION TECHNIQUE IN ORDER TO SOLVE THE RESULT
!  OF SUBSTITUTING EQUATIONS (4) AND (5) INTO (3),THEN (3) INTO (1) AND SOLVING.
!  FOR DOCUMENTATION, SEE REPORT "SOLVING FOR EEFF.SXW" UNDER REPORTS/MEMLS
!  FOLDER.
!
!  INPUTS: V - VOLUME FRACTION OF ICE
!          A - DEPOLARIZATION FACTOR
!          EICE - ICE PERMITTIVITY
!          EPSI - DRY SNOW PERMITTIVITY
!  OUTPUT: EPSEFF - EFFECITVE ICE PERMITTIVITY

IMPLICIT NONE

INTEGER, INTENT(IN) :: NUM
REAL, INTENT(IN) :: VFI(NUM),A(NUM),EICE,EPSI(NUM)
REAL, INTENT(OUT) :: EPSEFF(NUM)
REAL EPS,X,F,U1,V1,U2,V2,B,T,DBDX,DTDX,DFDX,DELTAX
INTEGER I,K

EPS=0.001
DO I=1,NUM
  X=1-EPSI(I)    !X0, INITIAL GUESS
  F=EPS*10       !DUMMY VALUE FOR FIRST GUESS
  K=1            !ITERATION COUNTER
  DO
    K=K+1
    U1=2*(1+(1-A(I))*X)
    V1=1+(1-A(I))*X+A(I)*(EICE-1)
    U2=1+2*A(I)*X
    V2=1+2*A(I)*X+(1-2*A(I))*(EICE-1)
    B=3-VFI(I)*(EICE-1)*(2*A(I)/V1+(1-2*A(I))/V2)
    T=VFI(I)*(EICE-1)*(U1/V1+U2/V2)
    DBDX=2*A(I)*VFI(I)*(EICE-1)*((1-A(I))/V1**2+(1-2*A(I))/V2**2)
    DTDX=DBDX*(EICE-1)
    DFDX=1-(B*DTDX-T*DBDX)/B**2
    DELTAX=-1*F/DFDX
    X=X+DELTAX
    F=X-T/B
    IF (ABS(F)<EPS) EXIT    ! TEST WHETHER WE ARE CLOSE ENOUGH TO ZERO
  END DO
  EPSEFF(I)=X+1
END DO
  
END SUBROUTINE POLDER

! -------------------------------------------------------------------------
!
SUBROUTINE POLMIX(TSCAT,SIH,SIV,NUM)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   CALCULATES THE POLARIZATION MIXING OF THE INTERFACE REFLECTIVITIES
!       OF EACH LAYER (TAKING INTO ACCOUNT THE FIRST ORDER SCATTERING)
!
!   [SIH,SIV] = POLMIX(TSCAT,SIH,SIV)
!       SIH:   INTERFACE REFLECTIVITY AT H-POL
!       SIV:   INTERFACE REFLECTIVITY AT V-POL
!       TSCAT: TAU SCAT
!
!   VERSION HISTORY:
!      1.0    WI 14.10.97
!      1.1    WI  4.11.97  BUG FIX (LAYER NUMBERING PROBLEM)
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!
!   USES - NONE
!
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER,INTENT(IN) :: NUM
REAL,INTENT(IN) :: TSCAT(NUM)
REAL,INTENT(INOUT) :: SIH(NUM+1),SIV(NUM+1)
REAL,DIMENSION(:),ALLOCATABLE ::  TSCATLONG,SMEAN,DELTAS

ALLOCATE( TSCATLONG(NUM+1),SMEAN(NUM+1),DELTAS(NUM+1) )

TSCATLONG(1:NUM)=TSCAT(1:NUM)
TSCATLONG(NUM+1)=1.

SMEAN=0.5*(SIH+SIV)
DELTAS=0.5*TSCATLONG*(SIH-SIV)

SIH=SMEAN+DELTAS
SIV=SMEAN-DELTAS

DEALLOCATE( TSCATLONG,SMEAN,DELTAS )

END SUBROUTINE POLMIX

! -------------------------------------------------------------------------
!
SUBROUTINE RO2EPSD(ROI,TI,FREQ,EPSI,EPSII,NUM,pixel,replicate,rank,meas,fn)
!
! -------------------------------------------------------------------------
!   
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   CALCULATES THE DIELECTRIC PERMITTIVITY FROM 
!   DENSITY FOR DRY SNOW.
!
!   [EPSI,EPSII] = RO2EPSD(ROI,TI,FREQ)
!       EPSI:  REAL PART OF DIELECTRIC PERMITTIVITY
!       EPSII: IMAGINARY PART OF DIELECTRIC PERMITTIVITY
!       ROI:   DENSITY
!       TI:    SNOW TEMPERATURE IN KELVIN
!       FREQ:  FREQUENCY
!
!   VERSION HISTORY:
!      1.0    WI 15.7.95
!      2.0    WI 12.11.97  ENHANCED WITH POLDER AND VAN SANTEN EQUATIONS (SEE
!                           POLDER.M)
!      3.0    MD  1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!      3.1    MD 21 NOV 05 CHANGED ALL LOCALS TO DYNAMICALLY ALLOCATABLE
!
!   USES:
!       EPSICE, EPSR, POLDER
!
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER,INTENT(IN):: NUM
REAL, INTENT(IN) :: ROI(NUM),TI(NUM),FREQ
REAL, INTENT(OUT) :: EPSI(NUM), EPSII(NUM)
INTEGER :: I
REAL :: EI
REAL,DIMENSION(:),ALLOCATABLE ::  EICE,F,A,EPSP,A3,EA,EA3,K1,K3,KSQ

integer,intent(in) :: pixel,replicate,rank,meas,fn

ALLOCATE( EICE(NUM),F(NUM),A(NUM),EPSP(NUM),A3(NUM),EA(NUM),EA3(NUM), & 
     K1(NUM),K3(NUM),KSQ(NUM) )

CALL EPSICE(TI,FREQ,EICE,NUM)

CALL EPSR(ROI,NUM,EPSI)

!     IMAGINARY PART AFTER POLDER AND VAN SANTEN 1946 (EFFECTIVE-MEDIUM APPROX)

F=ROI/0.917
EI=3.185

DO I=1,NUM
  A(I)=0.3
  IF (F(I)<0.55) THEN
    A(I)=0.476-0.64*F(I)
  END IF
  IF (F(I)<0.333) THEN
    A(I)=0.1+0.5*F(I)
  END IF
END DO

CALL POLDER(F,A,EI,EPSI,EPSP,NUM)

A3=1-2*A
EA=(EPSP*(1-A))+A
EA3=(EPSP*(1-A3))+A3
K1=(EA/(EA+A*(EI-1)))**2
K3=(EA3/(EA3+A3*(EI-1)))**2
KSQ=(2*K1+K3)/3
EPSII=EPSI**0.5*EICE*KSQ*F

DEALLOCATE( EICE,F,A,EPSP,A3,EA,EA3,K1,K3,KSQ )

END SUBROUTINE RO2EPSD

! -------------------------------------------------------------------------
!
SUBROUTINE RT(GAI,GBI,DEI,RI,TI,NUM)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   CALCULATES THE LAYER REFLECTIVITY AND TRANSMISSIVITY   
!
!
!   [RI,TI] = RT(GAI,GBI,DEI)
!       RI:   LAYER REFLECTIVITY
!       TI:   LAYER TRANSMISSIVITY
!       GAI:  ABSORPTION COEFFICIENT
!       GBI:  SCATTERING COEFFICIENT
!       DEI:  PATH LENGTH
!
!   VERSION HISTORY:
!      1.0    WI 15.7.95
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!   
!   USES: NONE
!
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER,INTENT(IN) :: NUM
REAL,INTENT(IN) :: GAI(NUM),GBI(NUM),DEI(NUM)
REAL,INTENT(OUT) :: RI(NUM),TI(NUM)
REAL,DIMENSION(:),ALLOCATABLE :: MYGAMMA,T0I,R0I,T02,R02
INTEGER I

ALLOCATE( MYGAMMA(NUM),T0I(NUM),R0I(NUM),T02(NUM),R02(NUM) )

MYGAMMA=(GAI*(GAI+2*GBI))**0.5
T0I=EXP(-1*MYGAMMA*DEI)

DO I=1,NUM
  IF(GBI(I)>0.00001)THEN
    R0I(I)=GBI(I)/(GAI(I)+GBI(I)+MYGAMMA(I))
  ELSE
    R0I(I)=0
  END IF
END DO

T02=T0I**2
R02=R0I**2
RI=R0I*(1-T02)/(1-T02*R02)
TI=T0I*(1-R02)/(1-T02*R02)

DEALLOCATE( MYGAMMA,T0I,R0I,T02,R02 )

END SUBROUTINE RT

! -------------------------------------------------------------------------
!
SUBROUTINE RUFFSOIL(F,MV,T,KSIGMA,THETA,EPS_TOP,&
             R_H_MOD,R_V_MOD)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM PULLIAINEN
!
!     FUNCTION FOR CALCULATING REFLECTIVITIES OF ROUGH, BARE SOILS
!     ACCORDING TO THEORY BY WEGMULLER & MÄTZLER
!     F = FREQUENCY [HZ]
!     MV = VOLUMETRIC MOISTURE [0..1]
!     T = TEMP [K]
!     KSIGMA = NORMALIZED SURFACE SDEV
!     THETA = NADIR ANGLE [DEG]
!     EPS_TOP = EPSILON OF OVERLYING MEDIUM
!     EPS_SOIL = DEFINE EPSILON EXPLICITLY
!
!    VERSION HISTORY:
!      1.0    ?? ?.?.?
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB


IMPLICIT NONE
REAL,INTENT(IN) :: F,MV,T,KSIGMA,THETA
COMPLEX, INTENT(IN) :: EPS_TOP
REAL,INTENT(OUT) :: R_H_MOD,R_V_MOD
REAL GND_TEMP,GND_EPS,PI,FRESNEL_H,FRESNEL_V,THETA_R
COMPLEX EPS_SOIL,EPS_EFF

PI=3.14159
GND_TEMP=T-273.15 ! K TO C

CALL EPSS(MV,GND_TEMP,F,EPS_SOIL)

EPS_EFF=EPS_SOIL/EPS_TOP

CALL GAMMAH(EPS_EFF,THETA,FRESNEL_H)
CALL GAMMAV(EPS_EFF,THETA,FRESNEL_V)

THETA_R=THETA/180*PI

R_H_MOD=FRESNEL_H*EXP(-KSIGMA**((0.1*COS(THETA_R))**0.5))

IF (THETA<=60) THEN
  R_V_MOD=R_H_MOD*COS(THETA_R)**0.65
ELSEIF(THETA==70)THEN
  R_V_MOD=R_H_MOD*0.621
ELSE
  PRINT*,("OOPS!  USE R_V_MOD=R_H_MOD, SINCE THETA=") ! ADDED BY MIKE
  PRINT *, THETA
  R_V_MOD=R_H_MOD 
END IF

END SUBROUTINE RUFFSOIL

! -------------------------------------------------------------------------
!
SUBROUTINE SCCOEFF(ROI,TI,PCI,FREQ,WIFR,GAI,SCCHO,GBIH,GBIV,GS6,&
                     GA2I,NUM,EPSI,EPSII,pixel,replicate,rank,meas)
!
! -------------------------------------------------------------------------
!
! CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER, MOD. BY MIKE (SEE HISTORY) 
!
!   CALCULATES THE SCATTERING COEFFICIENT FROM STRUCTURAL PARAMETERS
!     DIFFERENT ALGORITHMS CAN BE CHOSEN, BY CHANGING "SCCHO"
!
!   [GBIH,GBIV,GS6,GA2I] = SCCOEFF(ROI,TI,PCI,FREQ,WI,GAI,SCCHO)
!       GBIH:  2-FLUX SCATTERING COEFFICIENT AT H POL
!       GBIV:  2-FLUX SCATTERING COEFFICIENT AT V POL
!       GS6:   6-FLUX SCATTERING COEFFICIENT
!       GA2I:  2-FLUX ABSORPTION COEFFICIENT
!       ROI:   DENSITY
!       TI:    PHYSICAL TEMPERATURE
!       PCI:   CORRELATION LENGTH
!       FREQ:  FREQUENCY
!       WI:    WETNESS
!       GAI:   ABSORPTION COEFFICIENT
!       SCCHO: SCATTERING COEFFICIENT ALGORITHM CHOSEN
!
!   VERSION HISTORY:
!      1.0B    WI 15.7.95
!      1.0     WI 23.9.97 BUG FIXED
!      1.1     WI 26.9.97 LATEST FIT ON EXPERIMENTAL DATA WAS ADDED (OPTION 7)
!      1.2     WI 13.10.97 OPTION 8 ADDED, ADAPTED SCATTERING OF A SHELL/SPHERE
!                            TO NOTE 9/VER2 
!      1.3     WI  4.11.97 OPTION 9, 10 AND 11 ADDED 
!      1.4     WI 27.05.98 BORN APPROXIMATION ADDED (BORNA.M)
!      2.0     MD 1 APR 05 CALL TO BORNA CHANGED TO CALL TO BORN (BY MIKE). 
!                            ALL OPTIONS FOR OTHER SHAPES REMOVED BUT BORN 
!                            APPROX AND THE APPROX FOR FINE GRAIN SNOW IN 
!                            WIESMANN AND MATZLER 1999
!      2.0     MD 21 NOV 05  MADE ALL LOCALS DYNAMICALLY ALLOCATABLE
!
!   USES: BORN
!
!
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER, INTENT(IN) :: NUM,SCCHO
REAL, INTENT(IN) :: ROI(NUM),TI(NUM),PCI(NUM),FREQ,WIFR(NUM),&
                      GAI(NUM),EPSI(NUM),EPSII(NUM)
REAL, INTENT(OUT) :: GBIH(NUM),GBIV(NUM),GS6(NUM),GA2I(NUM)
REAL C,ROAIR,ROICE,DGB0H,DGB0V,K,EICE,PI
REAL,DIMENSION(:),ALLOCATABLE ::  VFI,KP,GB6,GC6,GF6,GTR,OMEGA
INTEGER :: I,IERR !IERR USED FOR MPI_FINALIZE IF WE HAVE AN ERROR

integer, intent(in) :: pixel,replicate,rank,meas

ALLOCATE( VFI(NUM),KP(NUM),GB6(NUM),GC6(NUM),GF6(NUM),GTR(NUM),OMEGA(NUM) )

! CONSTANTS
PI=3.14159
C=2.99
ROAIR=0.001293
ROICE=0.917

! SPECULAR COMPONENT OF SCATTERING COEFFICIENT
!   USUALLY 0 CAN BE IMPORTANT IN NEW SNOW!
DGB0H=0.
DGB0V=0.
K=FREQ*2.*PI/0.299793
EICE=3.18
VFI=ROI/ROICE

IF (SCCHO==1) THEN
  !DUE TO A BAD UDPATE, I'M ADDING A CHECK TO VERIFY THAT ROI IS REASONABLE
  !THE CONSTANTS IN THE FOLLOWING LINE ARE DERIVED FROM EMPIRICAL RELATIONSHIP
  !IN THE CALCULATION OF GS6
  DO I=1,NUM
    IF(ROI(I).GE.(7.48*PCI(I)+0.439))THEN
      PRINT *, 'ERROR IN SCCOEFF: DENSITY IS TOO LARGE OR GRAIN SIZE IS', &
        'TOO SMALL; ROI(I)=',ROI(I),'PCI(I)=',PCI(I),'LAYER#=',I,&
        '. STOPPING PROGRAM EXECUTION'
      CALL MPI_FINALIZE(IERR)
      STOP
    END IF
  END DO
  !FIT VOM 4.11.97 (WITH DENSITY, USES CORR. LENGTH FROM EXP. FIT!)
  GS6=(9.2*PCI-1.23*ROI+0.54)**2.5*(FREQ/50.)**2.5
  OMEGA=((EPSI-1)/EPSI)**0.5
  GB6=0.5*GS6*(1-OMEGA)
  GC6=0.25*GS6*OMEGA
ELSEIF (SCCHO==2) THEN
   ! BORN APPROXIMATION
   CALL BORN(K,VFI,PCI,EPSI,EICE,GB6,GC6,GF6,GS6,NUM,pixel,replicate,rank,meas)
END IF

! COMPUTE OTHER SCATTERING COEFFICIENTS

GTR=(4*GC6)/(GAI+2*GC6)
GA2I=GAI*(1+GTR)
GBIH=(GB6+DGB0H)+GTR*GC6
GBIV=(GB6+DGB0V)+GTR*GC6

DEALLOCATE( VFI,KP,GB6,GC6,GF6,GTR,OMEGA )

END SUBROUTINE SCCOEFF

! -------------------------------------------------------------------------
!
SUBROUTINE SLRED(NUM,ROI,EPSI,EPSII,TEI,SIH,SIV,DI,DEI,TI,PCI,WI,GAI,&
  FREQ,RNUM,XROI,XEPSI,XEPSII,XTEI,XSIH,XSIV,XDI,XDEI,XTI,XPCI,XWI,XGAI,&
  pixel,replicate,rank,meas)
!
! -------------------------------------------------------------------------
!
!     NOTE!! I DID NOT IMPLEMENT THE FULL SLRED SUBROUTINE FROM MEMLS
!       BECUASE I AM WORKING WITH THREE SNOW LAYERS.  THUS, IF THERE ARE
!       SUCEEDING COHERENT LAYERS WE KNOW THEY ARE LAYERS ONE AND TWO 
!       AND A HIGHLY SIMPLIFIED VERSION OF THE MATLAB CAN BE USED.

!       THIS CODE WILL NOT FUNCTION WITH SNOWPACKS WITH SUCCEEDING COHERENT 
!       LAYERS IF THERE ARE MORE THAN THREE LAYERS. OTHERWISE THE PROGRAM 
!       WILL STOP ALTOGETHER WITH AN ERROR MESSAGE. THIS CODE WILL CHECK WHETHER
!       THERE ARE COHERENT SCATTERING LAYERS AND INCORPORATE THESE INTO 
!       THE INTERFACE REFLECTIVITIES.  IF THERE ARE NO COHERENT LAYERS,
!       IT ASSIGNS THE INPUT ARRAYS TO THE OUTPUT ARRAYS.
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   LOCATES AND TREATS COHERENT LAYERS IN A SNOWPACK
!     SEE TECHNOTE 11
!
!   [RNUM,RROI,REPSI,REPSII,RTEI,RSIH,RSIV,RDI,RDEI,RTI,RPCI,RWI,RGAI]=
!       SLRED(NUM,ROI,EPSI,EPSII,TEI,SIH,SIV,DI,DEI,TI,PCI,FREQ,WI,GAI) 
!       NUM:  INDEX OF THE LAYER IN THE ORIGINAL SNOWPACK
!       ROI:  DENSITY [G/CM^3]
!       EPSI:
!       EPSII:
!       TEI:  LOCAL INCIDENCE ANGLE
!       SIH:  LAYER REFLECTIVITY AT H POL
!       SIV:  LAYER REFLECTIVITY AT V POL
!       DI:   LAYER THICKNESS
!       DEI:  LOCAL PATH LENGTH [M]
!       TI:   PHYSICAL SNOW TEMPERATURE [K]
!       PCI:  CORRELATION LENGTH [MM]
!       WI:   WETNESS  
!       GAI:  ABSORPTION COEFFICIENT
!       FREQ: FREQUENCY
!
!      1.0    WI 21.8.95
!      2.0    WI 13.8.98  COMPLETELY REWRITTEN
!      3.0    MD 1  APR 05 STRIPPED DOWN TO ONLY CHECK FOR COHERENT LAYERS -MD
!      3.1    MD 18 APR 05 FUNCTIONALITY ADDED TO DEAL WITH COHERENT LAYERS, 
!                            SO LONG AS THERE ARE NOT SUCCEEDING LAYERS
!      3.2    MD 04 NOV 05 FUNCTIONALITY ADDED TO DEAL WITH SUCCEEDING 
!                            COHERENT LAYERS FOR THREE LAYERS
!      3.3    MD 12 NOV 05 CHANGED TO DYNAMICALLY ALLOCATE ALL LOCALS
!   
!   USES:
!       FRESNELRC
!
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND

IMPLICIT NONE

INTEGER,INTENT(IN) :: NUM
REAL,INTENT(INOUT) :: ROI(NUM),EPSI(NUM),EPSII(NUM),TEI(NUM+1),SIH(NUM),&
  SIV(NUM),DI(NUM),DEI(NUM),TI(NUM),PCI(NUM),WI(NUM),GAI(NUM),FREQ
REAL,INTENT(OUT) :: XROI(NUM),XEPSI(NUM),XEPSII(NUM),XTEI(NUM+1),XSIH(NUM),&
  XSIV(NUM),XDI(NUM),XDEI(NUM),XTI(NUM),XPCI(NUM),XWI(NUM),XGAI(NUM)
INTEGER, INTENT(OUT) :: RNUM
INTEGER, DIMENSION(:), ALLOCATABLE :: I,A
REAL CC,PI,FIC,FC,THETA
REAL,DIMENSION(:),ALLOCATABLE ::  NS,FI,X
INTEGER J,NTHIN,PL,SC,SCMAX,ML,MLO,M,K,NAEQ0,NAEQ1,tnum
REAL FITOT
REAL,DIMENSION(:),ALLOCATABLE :: TROI,TEPSI,TEPSII,TTEI,FITV,&
  TSIH,TSIV,TDI,TDEI,TTI,TPCI,TWI,TGAI,TNS,TFI,EPSI_AUG,FH,FV

integer,intent(in) :: pixel,replicate,rank,meas

ALLOCATE( NS(NUM),FI(NUM), X(NUM) )

! IN GENERAL, WE FIRST TEST TO SEE IF THERE ARE ANY COHERENT LAYERS.  IF NOT,
!   WE SKIP TO CASE 2, OTHERWISE, WE HAVE CASE 1.  IN CASE 1, THE FIRST STEP IS
!   TO TREAT THE EFFECT OF SUCCEEDING COHERENT SNOWPACK LAYERS BY COMBINING 
!   THESE INTO A SINGLE LAYER. IF THERE ARE NO SUCCEEDING COHERENT LAYERS,
!   WE SKIP TO STEP 2; REGARDLESS, THE ARRAYS WITH THE T- PREFIX ARE DEFINED
!   BEFORE STEP 2.  IN STEP 2, THE EFFECT OF ISOLATED COHERENT LAYERS IS 
!   DEALT WITH BY INCORPORATING THE EFFECT OF EACH COHERENT LAYER INTO THE
!   INTER-LAYER REFLECTIVITIES.  

! CONSTANTS
CC=0.299793
PI=3.14159
FIC=4*PI*FREQ/CC
FC=4.712

! COMPUTE FI IN ORDER TO CHECK ON COHERENT LAYERS
THETA=TEI(NUM+1)
NS=EPSI**0.5
FI=FIC*DI*NS*COS(TEI(1:NUM))

! COUNT THE THIN LAYERS
NTHIN=0
DO J=1,NUM
  IF (FI(J)<FC) THEN
    NTHIN=NTHIN+1
  END IF
END DO

! ALLOCATE THE I ARRAY, WHICH CONTAINS THE INDECES OF THE THIN LAYERS
ALLOCATE(I(1:NTHIN),A(NUM))

! DEFINE THE I AND A ARRAYS.  THE A ARRAY HAS A VALUE OF 1 FOR THIN 
!   LAYERS AND 0 OTHERWISE
K=0
DO J=1,NUM
  IF (FI(J)<FC) THEN
    K=K+1
    I(K)=J
    A(J)=1
  ELSE
    A(J)=0
  END IF
END DO
! THE BOTTOM LAYER IS ALWAYS ASSUMED NONCOHERENT
A(1)=0

! CHECK TO SEE IF WE HAVE ANY COHERENT LAYERS TO TREAT.  IF NOT, ALLOCATE
!   THE REPLACEMENT ARRAYS (WITH THE 'R' PREFIX) WITH A SIZE OF THE ORIGINAL
!   NUMBER OF LAYERS, THEN DEFINE THOSE ARRAYS BASED ON THE ORIGINAL VALUES
! IF THERE IS ONLY ONE LAYER, IT IS ASSUMED TO BE NONCOHERENT, SO THE ORIGINAL
!   NUMBER OF LAYERS IS ASSIGNED TO RNUM AND THE OUTPUTS SET TO THE INPUTS.

IF (NTHIN>0.AND.NUM.GT.1) THEN
  ! CASE 1: TREAT COHERENT LAYERS IN THE SNOWPACK
  ! 1) IDENTIFY SUCCEEDING COHERENT LAYERS AND MARK THE PACKS FROM 2,3...SCMAX
  PL=0
  SC=1
  SCMAX=0
  ML=0
  MLO=0

  DO M=2,NUM
    IF (A(M)==1 .AND. PL==1) THEN
      IF(ML==0)THEN
        SC=SC+1
        ML=1
        A(M-1)=SC
      END IF

      A(M)=SC
      SCMAX=SC
    ELSE
      IF (PL==1) THEN
        PL=0
      ELSE
        IF(A(M)==1)THEN
          PL=1
          ML=0
        END IF
      END IF
    END IF
  END DO

  ! I AM NOT IMPLEMENTING THE GENERAL COMBINATION OF SUCCESSIVE COHERENT LAYERS.
  !   INSTEAD, I AM IMPLEMENTING A THREE-LAYER SOLUTION.  IF N~=3, THEN THE CODE
  !   WILL PRINT AN ERROR AND TERMINATE.
 
  IF (SCMAX>0) THEN
    IF(NUM.NE.3)THEN
      PRINT *,"THERE ARE SUCCESSIVE COHERENT LAYERS AND N~=3... ABORTING"
      STOP
    END IF
    ! THE FOLLOWING LINES ALLOCATE THE T- PREFIX ARRAYS BASED ON A TWO 
    ! LAYER PACK.  IF THIS POINT IN THE CODE IS REACHED, IT MEANS THAT BOTH
    ! BOTH LAYERS 2 AND 3 ARE COHERENT, SINCE LAYER 1 IS ASSUMED NON-COHERENT
    ! IN MEMLS IN GENERAL.  THEREFORE, A PHASE-WEIGHTED AVERAGE OF THE SNOW
    ! PROPERTIES IS COMPUTED AND THIS AVERAGE IS ASSIGNED TO LAYER 2 IN THE
    ! T- PREFIX ARRAY 

    FITOT=SUM(FI(2:3))
    ALLOCATE(FITV(2),TDI(2),TDEI(2),TROI(2),TTI(2),TWI(2),TPCI(2),TNS(2),&
      TGAI(2),TEPSI(2),TEPSII(2),TTEI(3),TFI(2))
    FITV=FI(2:3)/FITOT
    TDI(1)=DI(1)
    TDI(2)=SUM(DI(2:3))
    TDEI(1)=DEI(1)
    TDEI(2)=SUM(DEI(2:3))
    TROI(1)=ROI(1)
    TROI(2)=SUM(ROI(2:3)*FITV)
    TTI(1)=TI(1)
    TTI(2)=SUM(TI(2:3)*FITV)
    TWI(1)=WI(1)
    TWI(2)=SUM(WI(2:3)*FITV)
    TPCI(1)=PCI(1)
    TPCI(2)=SUM(PCI(2:3)*FITV)
    TNS(1)=NS(1)
    TNS(2)=SUM(NS(2:3)*FITV)
    TGAI(1)=GAI(1)
    TGAI(2)=SUM(GAI(2:3)*FITV)
    TEPSI(1)=EPSI(1)
    TEPSI(2)=SUM(EPSI(2:3)*FITV)
    TEPSII(1)=EPSII(1)
    TEPSII(2)=SUM(EPSII(2:3)*FITV)
    TTEI(1)=TEI(1)
    TTEI(2)=SUM(TEI(2:3)*FITV)
    TTEI(3)=TEI(4)
    TFI(1)=FI(1)
    TFI(2)=FITOT
    A(3)=1 !NOT SURE IF THIS JIVES WITH MY APPROACH NOW...
    TNUM=2
  ELSE
    !IF THERE ARE NO SUCCEEDING COHERENT LAYERS, I WILL MAP THE x ARRAYS TO 
    !  THE Tx ARRAYS...
    ALLOCATE(TDI(3),TDEI(3),TROI(3),TTI(3),TWI(3),TPCI(3),TNS(3),TGAI(3),&
      TEPSI(3),TEPSII(3),TTEI(4),TFI(3))
    TDI=DI
    TDEI=DEI
    TROI=ROI
    TTI=TI
    TWI=WI
    TPCI=PCI
    TNS=NS
    TGAI=GAI
    TEPSI=EPSI
    TEPSII=EPSII
    TTEI=TEI
    TFI=FI
    TNUM=3
  END IF
  DEALLOCATE(I,A)

  ! 2) DEAL WITH ISOLATED COHERENT LAYERS BY COMBINING THESE LAYERS' EFFECTS
  ! INTO THE INTER-LAYER REFLECTIVITIES
  
  ! COUNT THE THIN LAYERS AFTER REMOVAL OF SUCCEEDING COHERENT LAYERS
  NTHIN=0
  DO J=1,TNUM
    IF (TFI(J)<FC) THEN
      NTHIN=NTHIN+1
    END IF
  END DO

  ! ALLOCATE THE A ARRAY WHICH MARKS COHERENT LAYERS; THE I ARRAY SWITCHES
  !   BETWEEN MARKING INDECES OF COHERENT AND NONCOHERENT LAYERS.  THE GOAL
  !   IS TO CALCULATE THE INTERFACE REFLECTIVITIES.
  ALLOCATE(I(NTHIN),A(TNUM))

  ! REDEFINE THE I AND A ARRAYS.  THE A ARRAY HAS A VALUE OF 1 FOR THIN 
  !   LAYERS AND 0 OTHERWISE, NOW THAT THERE ARE NO SUCCEEDING COHERENT LAYERS
  K=0
  DO J=1,TNUM
    IF (TFI(J)<FC) THEN
      K=K+1
      I(K)=J
      A(J)=1
    ELSE
      A(J)=0
    END IF
  END DO
  ! THE BOTTOM LAYER IS ALWAYS ASSUMED NONCOHERENT
  A(1)=0
  DEALLOCATE(I) !I IS NOT USED NOW

  !  FIND THE VALUES OF A EQUAL TO ZERO, SAVE THE INDECES TO I
  NAEQ0=0
  DO J=1,TNUM
    IF (A(J)==0) THEN
      NAEQ0=NAEQ0+1
    END IF
  END DO
  ALLOCATE(I(NAEQ0))

  K=0
  DO J=1,TNUM
    IF (A(J)==0) THEN
      K=K+1
      I(K)=J
    END IF
  END DO

  ! CALCULATE INTERFACE REFLECTION COEFFICIENTS
  ALLOCATE(EPSI_AUG(TNUM+1),FH(TNUM),FV(TNUM))
  EPSI_AUG(1:TNUM)=TEPSI
  EPSI_AUG(TNUM+1)=1.0
  
  CALL FRESNELRC(TTEI,EPSI_AUG,FH,FV,TNUM)

  ! FOR LAYERS OF TYPE 0 SHI = FH^2
  ALLOCATE(TSIH(TNUM),TSIV(TNUM))
  DO J=1,NAEQ0
    TSIH(I(J))=FH(I(J))**2
    TSIV(I(J))=FV(I(J))**2
  END DO

  IF(NTHIN.EQ.0)THEN
    ! THERE ARE NO MORE COHERENT LAYERS AFTER REMOVAL OF SUCCEDING COHERENT 
    ! LAYERS SO SET THE OUTPUT ARRAYS EQUAL TO CURRENT Tx ARRAYS 
    RNUM=TNUM 
    XROI(1:RNUM)=TROI 
    XEPSI(1:RNUM)=TEPSI
    XEPSII(1:RNUM)=TEPSII
    XTEI(1:RNUM+1)=TTEI
    XSIH(1:RNUM)=TSIH
    XSIV(1:RNUM)=TSIV
    XDI(1:RNUM)=TDI
    XDEI(1:RNUM)=TDEI
    XTI(1:RNUM)=TTI
    XPCI(1:RNUM)=TPCI
    XWI(1:RNUM)=TWI
    XGAI(1:RNUM)=TGAI
  ELSE
    ! IF THERE ARE STILL COHERENT LAYERS AFTER REMOVAL OF SUCCEEDING COHERENT
    ! LAYERS, THEN DEAL WITH THESE.  DO A REDUCTION ON LAYERS OF TYPE 0 
    !   (COHERENT LAYER EFFECTS ARE TAKEN INTO ACCOUNT IN THE LAYER 
    !   REFLECTIVITIES)

    !  FIND THE NUMBER OF VALUES OF A EQUAL TO ONE

    DEALLOCATE(I)
    NAEQ1=0
    DO J=1,TNUM
      IF (A(J)==1) THEN
        NAEQ1=NAEQ1+1
      END IF
    END DO
    ALLOCATE(I(1:NAEQ1))

    !  FIND THE VALUES OF A EQUAL TO ONE
    K=0
    DO J=1,TNUM
      IF (A(J)==1) THEN
        K=K+1
        I(K)=J
      END IF
    END DO
  
   !COMPUTE NEW INTERFACE REFLECTIVITIES
    X=0 !INITIALIZE VARIABLES
!    XSIH=TSIH
!    XSIV=TSIV
    DO J=1,NAEQ1
      X(I(J))=2*FH(I(J))*FH(I(J)-1)*COS(FI(I(J)))
      TSIH(I(J)-1)=(FH(I(J))**2+FH(I(J)-1)**2+X(I(J)))/&
        (1+FH(I(J))**2*FH(I(J)-1)**2+X(I(J)))
      X(I(J))=2*FV(I(J))*FV(I(J)-1)*COS(FI(I(J)))
      TSIV(I(J)-1)=(FV(I(J))**2+FV(I(J)-1)**2+X(I(J)))/&
        (1+FV(I(J))**2*FV(I(J)-1)**2+X(I(J)))
    END DO

    !  FIND THE NUMBER OF VALUES OF A EQUAL TO ZERO
    DEALLOCATE(I)
    NAEQ0=0
    DO J=1,TNUM
      IF (A(J)==0) THEN
        NAEQ0=NAEQ0+1
      END IF
    END DO
    ALLOCATE(I(NAEQ0))

    !  FIND THE VALUES OF A EQUAL TO ZERO
    K=0
    DO J=1,TNUM
      IF (A(J)==0) THEN
        K=K+1
        I(K)=J
      END IF
    END DO

    !  DEFINE OUTPUT VARIABLES BASED ON I ARRAY
    DO J=1,NAEQ0
!      XROI(J)=ROI(I(J))
!      XTEI(J)=TEI(I(J))
!      XDI(J)=DI(I(J))
!      XTI(J)=TI(I(J))
!      XPCI(J)=PCI(I(J))
!      XWI(J)=WI(I(J))
!      XGAI(J)=GAI(I(J))
      XROI(J)=TROI(I(J))
      XTEI(J)=TTEI(I(J))
      XDI(J)=TDI(I(J))
      XTI(J)=TTI(I(J))
      XPCI(J)=TPCI(I(J))
      XWI(J)=TWI(I(J))
      XGAI(J)=TGAI(I(J))
   
!      XSIH(J)=XSIH(I(J))
!      XSIV(J)=XSIV(I(J))

      XSIH(J)=TSIH(I(J))
      XSIV(J)=TSIV(I(J))

!      XEPSI(J)=EPSI(I(J))
!      XEPSII(J)=EPSII(I(J))
!      XDEI(J)=DEI(I(J))
      XEPSI(J)=TEPSI(I(J))
      XEPSII(J)=TEPSII(I(J))
      XDEI(J)=TDEI(I(J))
    END DO

    RNUM=TNUM-NTHIN
    XTEI(RNUM+1)=THETA 
  END IF
  DEALLOCATE(TDI,TDEI,TROI,TTI,TWI,TPCI,TNS,TGAI,TEPSI,TEPSII,TTEI,TFI,&
     TSIH,TSIV)
ELSE
  ! THERE ARE NO COHERENT LAYERS (OR ONLY ONE LAYER), SO SET THE OUTPUT 
  ! ARRAYS EQUAL TO THE INPUT ARRAYS
  RNUM=NUM 
  XROI=ROI
  XEPSI=EPSI
  XEPSII=EPSII
  XTEI=TEI
  XSIH=SIH
  XSIV=SIV
  XDI=DI
  XDEI=DEI
  XTI=TI
  XPCI=PCI
  XWI=WI
  XGAI=GAI
END IF

DEALLOCATE( NS,FI, X )
 
END SUBROUTINE SLRED

! -------------------------------------------------------------------------
!
SUBROUTINE SNOWAO(V,A,NUM)
!
! -------------------------------------------------------------------------
!
!     CODE ORIGINALLY OBTAINED IN MATLAB FROM MATZLER
!
!   COMPUTES THE DEPOLARIZATION FACTOR OF OBLATE SNOW GRAINS
!   NOTE 10, MÄTZLER 1997
!
!   A = SNOWAO(V)
!       A:    DEPOLARIZATION FACTOR OF OBLATE SNOW GRAINS
!       V:    VOLUME FRACTION OF ICE
!
!   VERSION HISTORY:
!      1.0    WI 29.5.98
!      2.0    MD 1 APR 05 TRANSLATED TO FORTRAN FROM MATLAB
!
!   USES: NONE
!       
!   COPYRIGHT (C) 1997 BY THE INSTITUTE OF APPLIED PHYSICS, 
!   UNIVERSITY OF BERN, SWITZERLAND


IMPLICIT NONE
INTEGER, INTENT(IN) :: NUM
REAL, INTENT(IN) :: V(NUM)
REAL, INTENT(OUT) :: A(NUM)
INTEGER I

DO I=1,NUM
  A(I)=0.3
  IF (V(I)<0.55) THEN
    A(I)=0.476-0.64*V(I)
  END IF
  IF (V(I)<=0.333) THEN
    A(I)=0.1+0.5*V(I)
  END IF
END DO

END SUBROUTINE SNOWAO