interfaez.pixel

! ----------------------------------------------------------------------
!
SUBROUTINE INTERFACEZ(Y,ZOUT,N_RL,N_R,N_Ypp,N_rP,N_Yr,N_Z,N_ZP,N_ZR,N_C,N_X,&
  CTRL,FREQ,THETA,X,N_A,N_U,FORCING,MONTH,TBRAW,ALBEDO,MEAS,RANK,VEG_MAP,&
  IERR,F_VEG,BDRF,N_BDRF,MEAS_SWITCH,VCOVER)
!
! ----------------------------------------------------------------------
!
! THIS SUBROUTINE PERFORMS SEVERAL FUNCTIONS.  OVER THE LOOP OF ALL THE 
! REPLICATES, 1) THE PREDICTED OBSERVATIONS ARE CALCULATED AT THE STATE
! PIXEL SCALE 2) THE PREDICTED OBSERVATIONS ARE AGGREGATED TO THE 
! MEASUREMENT SCALE FOR COMPARISON WITH THE ACTUAL OBSERVATIONS
!
! Y - STATE VECTOR
! Z - OUTPUT PREDICTED OBSERVATIONS
! N_RL - NUMBER OF REPLICATES LOCAL ON THE CALLING PROCESS; THIS ARGUMENT IS
!   SET TO 1 WHEN CALLING FROM ENKFGEN, SINCE WE ARE ONLY COMPUTING A SINGLE 
!   PREDICTED OBSERVATION
! N_R - NUMBER OF REPLICATES GLOBALLY
! N_Y - TOTAL NUMBER OF STATES
! N_rP - TOTAL NUMBER OF running PIXELS AT STATE RESOLUTION
! N_Z - TOTAL NUMBER OF OBSERVATIONS
! N_ZP - ARRAY OF TOTAL NUMBER OF MEASUREMENTS FOR EACH RESOLUTION
! N_ZR - NUMBER OF MEASUREMENTS (AT DIFFERENT RESOLUTIONS) OF EACH PIXEL 
! N_C - LENGTH OF CONTROL ARRAY (CTRL)
! N_X - DIMENSION OF X VARIABLE
! N_A - DIMENSION OF ALPHA VARIABLE
! CTRL - THE ARRAY WITH PARAMETERS CONTROLLING SIMULATIONS
! FREQ - FREQUENCY VECTOR
! THETA - MEASUREMENT ANGLE VECTOR
! X - AUXILIARY STATE VARIABLES
! N_U - DIMENSION OF FORCING VARIABLE 
! FORCING - FORCING DATA FOR STATE MODEL, NEEDED FOR RTMs
! MONTH AT WHICH RTM IS RUN
! TBRAW - PIXEL BRIGHTNESS TEMPS TO BE AGGREGATED
! ALBEDO - RAW ALBEDO OUTPUT
! MEAS - MEASUREMENT NUMBER
! RANK - MPI PROCESS RANK
! VEG_MAP - THE SSIB VEGETATION COVER CLASSIFICATION FOR EACH PIXEL
! IERR - MPI ERROR VARIABLE
! F_VEG - FACTOR BY WHICH VEGETATION IS SCALED (SPATIALLY UNIFORM)
! BDRF - ARRAY CONTAINING THE ALBEDO LOOKUP TABLE
! N_BDRF - NUMBER OF ROWS (ENTRIES) IN BDRF LUT
! MEAS_SWITCH - SWITCH THAT SPECIFIES WHICH MEASUREMENTS ARE USED. IF EQUAL
!   TO 3, THEN NO PM MEASUREMENTS ARE CALCULATED
! VCOVER - VEGETATION COVER PERCENTAGE
!
! CALLS: SSIB_LAYER, RAD_XFER_MOD
! CALLED FROM: ENKF, ENKFGEN
! 
! VERSION HISTORY: WRITTEN - MD - MAY 2005
!                  SET UP FOR VEG UNCERTAINTY - MD - JUN 2006

IMPLICIT NONE

INTEGER,INTENT(IN)::N_Ypp,N_rP,N_ZR,N_RL,N_R,N_Z,N_C,N_X,N_A,N_U,MONTH,N_YR,RANK,&
  IERR,N_BDRF,MEAS_SWITCH
INTEGER,INTENT(IN)::N_ZP(N_ZR),VEG_MAP(N_rP)
INTEGER,INTENT(INOUT) :: CTRL(N_C)
REAL,INTENT(IN)::Y(N_Ypp,N_RL),FREQ(CTRL(9)),THETA(CTRL(9)),X(N_X*N_rP,N_RL),&
  FORCING(N_U,N_RL,N_rP),ALBEDO(N_rP,N_RL), F_VEG(N_R,1),BDRF(N_BDRF,2),&
  VCOVER(N_rP)
REAL,INTENT(INOUT)::ZOUT(1:N_Z,N_RL),TBRAW(N_rP*N_ZP(1),N_RL)
INTEGER R,P,YI(3),ZI(4),i,J,K,L,M,N,O,S,T,KGLOBAL
REAL, DIMENSION(:), ALLOCATABLE :: DZ,ATM_IN,CAN_IN,AUX_SNOW_IN
REAL,DIMENSION(:,:),ALLOCATABLE ::  SNOW_IN,ZHAT,TB
INTEGER,INTENT(IN) :: MEAS
INTEGER,PARAMETER :: N_VEG=13 !NUMBER OF SSIB VEGETATION TYPES
INTEGER,PARAMETER :: N_MOS=12 !NUMBER OF MONTHS

!SSIB AND SLA VEGETATION DATA
REAL :: POROSITY(N_VEG), LAI(N_VEG,N_MOS,2),SLA(N_VEG),TRAN_VEG_NIR(N_VEG),&
  REF_VEG_NIR(N_VEG)
REAL :: RHO_VEG !USED FOR VEGETATION INPUTS

REAL :: RADIUS_IN,R_SNOW !USED FOR INTERPOLATING ALBEDOS AND INTERPOLATED OUTPUT
REAL :: R_VEG,T_VEG,V_C !USED IN NIR RTM

ALLOCATE( ZHAT(N_rP,2),DZ(3),ATM_IN(CTRL(8)),CAN_IN(CTRL(7)),&
  AUX_SNOW_IN(CTRL(5)), TB(2,CTRL(9)) )

! DEFINE CTRL PARAMETERS
CTRL(1)=N_rP !TOTAL NUMBER OF CALCULATIONS TO DO IS NUMBER OF STATE PIXELS
!NOTE: ATM RTM SWITCH CTRL(3) NOW SET IN RUN_PARAMS
!NOTE: CANOPY RTM SWITCH CTRL(4) NOW VARIES WITH DATA

! CALL TO GET SSIB VEGETATION DATA 
CALL GET_VEG_DATA(POROSITY,LAI,TRAN_VEG_NIR,REF_VEG_NIR)

! CALL TO GET SPECIFIC LEAF AREA DATA
CALL GET_SLA_DATA(SLA)
! DEFINE VEGETATION DENSITY
RHO_VEG=950.

DO R=1,N_RL

  CTRL(2)=3!OBTAIN NUMBER OF SNOW LAYERS FROM AUXILIARY SNOW ARRAY  
  ZOUT(1:12,R)=0. !INITIALIZE ZOUT SUMS FOR COMPUTING AVERAGE TB
  IF(N_R.EQ.1)THEN
    KGLOBAL=1 !IF THERE IS ONE GLOBAL REPLICATE (I.E., GENERATOR IS RUNNING) 
  ELSE
    KGLOBAL=RANK*N_RL+R !GLOBAL REPLICATE NUMBER, FOR INDEXING F_VEG FOR FILTER
  END IF

  DO P=1,N_rP

    ! 0) ISOLATE VEGETATION COVER FRACTION FOR THIS PIXEL
    V_C=VCOVER(P)/100.

    ! 1) SNOW AND SOIL INPUTS
    J=N_YR*(P-1)+1 !BEGINNING Y INDEX
    K=N_YR*P       !ENDING Y INDEX
    L=N_X*(P-1)+1  !BEGINNING X INDEX
    M=N_X*P        !ENDING X INDEX
    N=N_ZP(1)*(P-1)+1 !BEGINNING TB RAW INDEX
    O=N_ZP(1)*P       !ENDING TB RAW INDEX
    ! DETERMINE NUMBER OF SNOW LAYERS AND ASSIGN SNOW_IN ARRAY

    IF(Y(J,R).LT.1e-4)THEN !IF THERE IS LESS THAN 0.1MM, ASSUME NO SNOW
      CTRL(2)=1 !NO SNOW; THIS IS SET TO 1 FOR ALLOCATING SNOW_IN
      AUX_SNOW_IN(1)=0 ! NO SNOW
      ALLOCATE(SNOW_IN(CTRL(2),CTRL(6)))
      SNOW_IN=0. !FOR NO SNOW CASE, SET THE SNOW_IN ARRAY TO A DUMMY VALUE 0.
    ELSEIF(Y(J,R).LT.0.05)THEN
      CTRL(2)=1 ! ONE LAYER
      AUX_SNOW_IN(1)=1 ! ONE LAYER
      ALLOCATE(SNOW_IN(CTRL(2),CTRL(6)))
      SNOW_IN(1,1)=Y(J,R) !SNOW DEPTH
      SNOW_IN(1,2)=Y(J+1,R)  !ONE LAYER SNOW DENSITY
      SNOW_IN(1,3)=Y(J+10,R) !ONE LAYER SNOW GRAIN SIZE
      SNOW_IN(1,4)=Y(J+8,R)  !ONE LAYER SNOW LIQUID WATER CONTENT
      SNOW_IN(1,5)=Y(J+5,R)  !ONE LAYER SNOW TEMPERATURE
    ELSE
      CTRL(2)=3 !THREE SNOW LAYERS
      AUX_SNOW_IN(1)=3 ! THREE SNOW LAYERS
      ALLOCATE(SNOW_IN(CTRL(2),CTRL(6)))
      CALL SSIB_LAYER(Y(J,R),DZ) 
      SNOW_IN(1:3,1)=DZ(1:3)        !THREE LAYER SNOW DEPTH
      SNOW_IN(1:3,2)=Y(J+1:J+3,R)   !THREE LAYER SNOW DENSITY
      SNOW_IN(1:3,3)=Y(J+10:J+12,R) !THREE LAYER SNOW GRAIN SIZE
      SNOW_IN(1:3,4)=Y(J+7:J+9,R)   !THREE LAYER SNOW LIQUID WATER CONTENT
      SNOW_IN(1:3,5)=Y(J+4:J+6,R)   !THREE LAYER SNOW TEMPERATURE

!      if(meas.eq.85) then
!        print *, 'snow profile for m=85'
!        print *, snow_in(1,1:5)
!        print *, snow_in(2,1:5)
!        print *, snow_in(3,1:5)
!      end if

    END IF

    AUX_SNOW_IN(2)=Y(J+13,R)           !SOIL TEMPERATURE
    AUX_SNOW_IN(3)=X(L+4,R)            !SOIL MOISTURE (USE UPPER LAYER)
    AUX_SNOW_IN(4)=POROSITY(VEG_MAP(P))!SOIL POROSITY
    AUX_SNOW_IN(5)=0.16                !CONSTANT OF PROPORTIONALITY BETWEEN
                                       !GRAIN SIZE AND CORRELATION LENGTH

    ! 2) DEFINE CANOPY PM RTM PROPERTIES
    IF(LAI(VEG_MAP(P),MONTH,1).LT.1.1E-3.OR.SLA(VEG_MAP(P)).EQ.0.0 &
         .OR.SLA(VEG_MAP(P)).EQ.-9999.) THEN
      ! IF LEAF AREA INDEX IS LESS THAN 0.001, OR SLA IS 0.0, THEN DON'T
      !   APPLY THE VEGETATION MODEL
      CTRL(4)=0
    ELSE
      CTRL(4)=1                 !TURN ON VEGETATION RTM SWITCH
      CAN_IN(1)=8               !VEGETATION WATER SALINITY: [PPT]
      CAN_IN(2)=X(L+1,R)          !CANOPY TEMPERATURE [k]
      CAN_IN(3)=0.5             !VEGETATION GRAVIMETRIC WATER CONTENT [FRAC]
      ! NEEDLE DIAMETER (M) IS CALCULATED FROM SPECIFIC LEAF AREA AND VEGETATION
      !   DENSITY
      CAN_IN(4)=1./SLA(VEG_MAP(P))/RHO_VEG
      ! VEGETATION BIOMASS (KG VEGETATION /M2 GROUND AREA) IS CALCULATED AS THE 
      !   QUOTIENT OF LEAF AREA INDEX ( M2 VEGETATION AREA  / M2 GROUND AREA) 
      !   AND SPECIFIC LEAF AREA (M2 VEGETATION AREA / KG VEGETATION MASS)
      ! LEAF AREA INDEX PERTURBATION (F_VEG) FACTOR IS APPLIED HERE.
      CAN_IN(5)=LAI(VEG_MAP(P),MONTH,1)*EXP(F_VEG(KGLOBAL,1))/SLA(VEG_MAP(P))
    END IF
  
!    print *, ctrl(4)

    ! 3) DEFINE ATMOSPHERIC RTM PROPERTIES
    ATM_IN(1)=FORCING(5,R,P)/100.!CONVERT SURFACE PRESSURE TO MBAR FROM pA
    ATM_IN(2)=FORCING(6,R,P)     !SURFACE ATMOSPHERIC TEMPERATURE 
    ATM_IN(3)=MONTH              !MONTH OF YEAR
    ATM_IN(4)=1200*FORCING(7,R,P)!CONVERT MOISTURE TO G/M3 FROM G/G : ASSUME 
                                 ! AN AIR DENSITY OF 1200 G/M3

    ! 4) COMPUTE PREDICTED BRIGHTNESS TEMPERATURES AT EACH PIXEL ONLY 
    !      IF MEASUREMENT INTERVAL IS ODD, WHICH CORRESPONDS TO THE 1 am
    !      MEASUREMENTS. DO NOT COMPUTE IF THE SWITCH INDICATES NIR 
    !      MEASUREMENTS ONLY, OR IF IT INDICATES NO MEASUREMENTS
 
! THIS OPTION OUGHT TO BE AUTOMATED: DURING THE SYNTHETIC TESTS, I ALWAYS
!   CALCULATED ONLY THE NIGHTTIME PM MEASUREMENTS, WHICH SAVED ALOT OF TIME.
!   THIS OPTION CAN BE DISABLED BY RE-ARRANGING THE COMMENTS IN THE 'IF' BLOCK
!   BELOW, WHICH I DID FOR THE CLPX GBMR WORK.

!    IF(MOD(MEAS,2).EQ.0)THEN
!      ! DON'T COMPUTE BRIGHTNESS TEMPS DURING THE DAYTIME MEASUREMENTS
!      TB=-9999.
!    ELSEIF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
    IF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
!    IF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
      ! THEN WE ARE ONLY USING ALBEDO MEASUREMENTS OR NOT UPDATING THIS TIME
      TB=-9999
    ELSE
      CALL RAD_XFER_MOD(CTRL,FREQ,THETA,AUX_SNOW_IN,SNOW_IN,CAN_IN,ATM_IN,&
        TB,P,R,RANK,MEAS,V_C,IERR) 
    END IF

    ! 5) AGGREGATE PREDICTED OBSERVATIONS OVER THE DOMAIN
    !AGGREGATE MEASUREMENT 1 TO MEASUREMENT SCALE

    !RECORD TB RAW: SEE JOURNAL ON 10/10/2005 ON INDEXING
    DO I=1,CTRL(9)
      TBRAW(N+I-1,R)=TB(1,I)
      TBRAW(CTRL(9)+N+I-1,R)=TB(2,I)
    END DO

! I AM COMMENTING OUT THE PART THAT ONLY CALCULATES TB FOR DAYTIME MEAS...
    !ADD HORIZONTAL TB TO ZOUT
    DO I=1,CTRL(9)
!      IF(MOD(MEAS,2).EQ.0)THEN
!        !IT'S DAYTIME
!        ZOUT(I,R)=-9999.
!      ELSEIF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
      IF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
        ! THEN WE ARE ONLY USING ALBEDO MEASUREMENTS OR NOT UPDATING THIS TIME
        ZOUT(I,R)=-9999.
      ELSE
        ZOUT(I,R)=ZOUT(I,R)+TB(1,I)/N_rP !THIS WILL BE AN AVERAGE
      END If
    END DO
    !ADD VERTICAL TB TO ZOUT
    DO I=1,CTRL(9)
!      IF(MOD(MEAS,2).EQ.0)THEN
!        !IT'S DAYTIME
!        ZOUT(I+6,R)=-9999.
!      ELSEIF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
      IF(MEAS_SWITCH.EQ.3.OR.MEAS_SWITCH.EQ.9)THEN
        ! THEN WE ARE ONLY USING ALBEDO MEASUREMENTS OR NOT UPDATING THIS TIME
        ZOUT(I+6,R)=-9999.
      ELSE
        ZOUT(I+6,R)=ZOUT(I+6,R)+TB(2,I)/N_rP !THIS WILL BE AN AVERAGE
      END IF
    END DO

    !ALBEDOS...
!    THIS WAS THE OLD BROADBAND MEASUREMENT
!    ZOUT(N_ZP(1)+P,R)=ALBEDO(P,R) !ALBEDO FROM SSIB...

    !6) COMPUTE NARROWBAND ALBEDOS
    IF(MOD(MEAS,2).EQ.0)THEN
      !IT'S DAYTTIME

      !ISOLATE VEGETATION NIR RTM PROPERTIES FOR THIS PIXEL
      T_VEG=TRAN_VEG_NIR(VEG_MAP(P))
      R_VEG=REF_VEG_NIR(VEG_MAP(P))

      IF ( CTRL(2).EQ.3 )THEN
        RADIUS_IN=0.5*SNOW_IN(3,3)*1.0E6
        CALL INTERPOLATE(BDRF(1:N_BDRF,1),BDRF(1:N_BDRF,2),RADIUS_IN,&
          R_SNOW,N_BDRF)
          !FRACTIONAL VEGETATION COVERAGE MEASUREMENT RELATIONSHIPS 
          ZOUT(N_ZP(1)+P,R)=R_SNOW*(1-V_C)+(R_SNOW*T_VEG**2+R_VEG)*V_C
      ELSEIF ( CTRL(2).EQ.1 )THEN
        IF( AUX_SNOW_IN(1).EQ.1) THEN
          RADIUS_IN=0.5*SNOW_IN(1,3)*1.0E6
          CALL INTERPOLATE(BDRF(1:N_BDRF,1),BDRF(1:N_BDRF,2),RADIUS_IN,&
            R_SNOW,N_BDRF)
          !FRACTIONAL VEGETATION COVERAGE MEASUREMENT RELATIONSHIPS 
          ZOUT(N_ZP(1)+P,R)=R_SNOW*(1-V_C)+(R_SNOW*T_VEG**2+R_VEG)*V_C
        ELSEIF (AUX_SNOW_IN(1).EQ.0) THEN
          !if there is no snow, set albedo to 0.1?
          ZOUT(N_ZP(1)+P,R)=0.1
        ELSE
          PRINT *, 'INTERFACEZ: ERROR IN ALBEDO ASSIGNMENT CASE STRUCTURE 1.',&
            'ABORTING'
          CALL MPI_FINALIZE(IERR) 
          STOP
        END IF
      ELSE 
        PRINT *, 'INTERFACEZ: ERROR IN ALBEDO ASSIGNMENT CASE STRUCTURE 2.',&
          'ABORTING'
        CALL MPI_FINALIZE(IERR) 
        STOP
      END IF
    ELSE 
      !IT'S NIGHTTIME, SET BDRF TO NODATA 
      ZOUT(N_ZP(1)+P,R)=-9999.
    END IF


    !7) SURFACE TEMPERATURE MEASUREMENT
    IF(CTRL(4).EQ.1) THEN
      !OLD WAY:  ZOUT(SUM(N_ZP(1:2))+P,R)=X(L+1,R) 

      !IF THE VEGETATION RTM SWITCH IS ON CTRL(4)=1), THEN WE MUST TAKE THE 
      !  VEGETATION TEMPERATURE INTO ACCOUNT.  APPLY SIMPLE FRACTIONAL 
      !  FRACTIONAL VEGETATION COVERAGE MODEL

      IF(Y(J,R).LT.1e-4)THEN !IF THERE IS LESS THAN 0.1MM, ASSUME NO SNOW
        !USE GROUND TEMPERATURE
        ZOUT(SUM(N_ZP(1:2))+P,R)=V_C*X(L+1,R)+(1.-V_C)*AUX_SNOW_IN(2)
      ELSEIF(Y(J,R).LT.0.05)THEN
        !USE ONE LAYER SNOW TEMPERATURE
        ZOUT(SUM(N_ZP(1:2))+P,R)=V_C*X(L+1,R)+(1.-V_C)*SNOW_IN(1,5)
      ELSE
        !USE SURFACE SNOW TEMPERATURE
        ZOUT(SUM(N_ZP(1:2))+P,R)=V_C*X(L+1,R)+(1.-V_C)*SNOW_IN(3,5)
      END IF
    ELSEIF(CTRL(4).EQ.0) THEN
      !THIS CODE IS NOW REFERENCE BASICALLY IF WE HAVE BARE SOIL
      IF(Y(J,R).LT.1e-4)THEN !IF THERE IS LESS THAN 0.1MM, ASSUME NO SNOW
        !USE GROUND TEMPERATURE
        ZOUT(SUM(N_ZP(1:2))+P,R)= AUX_SNOW_IN(2)
      ELSEIF(Y(J,R).LT.0.05)THEN
        !USE ONE LAYER SNOW TEMPERATURE
        ZOUT(SUM(N_ZP(1:2))+P,R)= SNOW_IN(1,5)
      ELSE
        !USE SURFACE SNOW TEMPERATURE
        ZOUT(SUM(N_ZP(1:2))+P,R)= SNOW_IN(3,5)
      END IF
    ELSE
      PRINT *, 'ERROR IN INTERFACEZ; ILLEGAL VALUE OF CTRL(4). STOPPING.'
      CALL MPI_FINALIZE(IERR)
      STOP
    END IF

    DEALLOCATE(SNOW_IN)
  END DO
END DO

deallocate( ZHAT,dz,atm_in,can_in,aux_snow_in,tb)

END SUBROUTINE INTERFACEZ