!+ Program to set aerosol mixing ratios for SRA profiles.
!
! Method:
!	The pressure, temperatures and specific humidity are read
!	in as SAMSON files. These quantities enable the height at
!	each pressure level to be calculated: the SRA profiles are
!	specified in terms of heights. For the profile selected
!	the extinction at 550 nm for each component is calculated
!	and the volume fraction is set to be consistent with the
!	extinction.
!
! Current Owner of Code: J. M. Edwards
!
! History:
!	Version		Date		Comment
!	1.0		12-04-95	First Version under RCS
!					(J. M. Edwards)
!
! Description of Code:
!   FORTRAN 77 with extensions listed in documentation.
!
!- ---------------------------------------------------------------------
      PROGRAM SRA_AEROSOL_PROFILE
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_spec.ucf'
      INCLUDE 'incl/dimensions_field.ucf'
      INCLUDE 'incl/physical_constants_1.ccf'
      INCLUDE 'incl/gas_list.pcf'
      INCLUDE 'incl/aerosol_component.pcf'
      INCLUDE 'incl/aerosol_component.ncf'
      INCLUDE 'incl/aerosol_profile.ncf'
      INCLUDE 'incl/aerosol_profile.pcf'
      INCLUDE 'incl/aerosol_model.ncf'
      INCLUDE 'incl/aerosol_model.pcf'
      INCLUDE 'incl/aerosol_representation.pcf'
      INCLUDE 'incl/input_head.pcf'
      INCLUDE 'incl/input_head.ncf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/prec_machine.ucf'
      INCLUDE 'incl/precision.pcf'
      INCLUDE 'incl/error.pcf'
!
!     DECLARATION OF VARIABLES.
      CHARACTER
     &     BASE_NAME*80
!             BASE NAMES OF INPUT FILES
     &   , FILE_NAME*80
!             NAME OF CURRENT FILE
     &   , TITLE_DUMMY*80
!             DUMMY TITLE
      INTEGER
     &     IERR
!             ERROR FLAG
     &   , IOS
!             I/O ERROR FLAG
!
      LOGICAL
     &     LOCK_CODE
!             LOGICAL TO FORBID INTERACTIVE LOOPING
      EXTERNAL
     &     LOCK_CODE
!
      INTEGER
     &     N_LAYER
!             NUMBER OF LAYERS
     &   , N_LAYER_DUMMY
!             DUMMY VARIABLE
     &   , NDUM
!             DUMMY VARIABLE
     &   , I_AEROSOL_PROFILE
!             INDEX OF AEROSOL PROFILE
      REAL*8
     &     Z(0: NPD_LAYER)
!             HEIGHTS OF LEVELS
     &   , P(0: NPD_LAYER)
!             PRESSURES
     &   , T(0: NPD_LAYER)
!             TEMPERATURES
     &   , Q(0: NPD_LAYER)
!             SPECIFIC HUMIDITY
     &   , DENSITY(0: NPD_LAYER)
!             DENSITIES AT LEVELS
     &   , AEROSOL_MIX_RATIO(0: NPD_LAYER
     &   , NPD_AEROSOL_SPECIES)
!             NUMBER DENSITY OF AEROSOL
     &   , PSTAR
!             SURFACE PRESSURE
     &   , TSTAR
!             SURFACE AIR TEMPERATURE
     &   , QSTAR
!             SURFACE HUMIDITY
     &   , DENSITY_STAR
!             SURFACE DENSITY
     &   , P_DUMMY(0: NPD_LAYER)
!             DUMMY PRESSURE ARRAY
     &   , PSTAR_DUMMY
!             DUMMY SURFACE PRESSURE
      INTEGER
     &     TYPE_AEROSOL(NPD_AEROSOL_COMPONENT)
!             TYPES OF AEROSOLS
     &   , I
!             LOOP INDEX
     &   , J
!             LOOP INDEX
     &   , I_MODEL
!             INDEX OF AEROSOL MODEL
     &   , I_AEROSOL_LAYER
!             INDEX OF AEROSOL LAYER
     &   , LENGTH_NAME
!             LENGTH OF BASE NAME
     &   , I_COMPONENT
!             TEMPORARY INDEX
      REAL*8
     &     VOLUME_FRACTION(0: NPD_LAYER, NPD_AEROSOL_COMPONENT)
!             CALCULATED AEROSOL FRACTIONS
     &   , VOLUME_FRACTION_LEVEL(NPD_AEROSOL_COMPONENT)
!             FRACTIONS OF COMPONENTS AT LEVELS
     &   , EXT_550NM_OVERALL
!             OVERALL EXTINCTION AT 550NM
!
!     FUNCTIONS CALLED:
      REAL*8
     &     EXTINCTION_PROFILE
!             FUNCTION TO CALCULATE EXTINCTION
     &   , FNC_DENSITY
!             FUNCTION TO CALCULATE DENSITY
      EXTERNAL
     &     EXTINCTION_PROFILE, FNC_DENSITY
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     READ_SAMSON_FIELD, CALC_VOLUME_FRACTION
     &   , WRITE_SAMSON_FIELD
!
!
      INCLUDE 'incl/initial.pcf'
      INCLUDE 'incl/aerosol_profile.dcf'
      INCLUDE 'incl/aerosol_model.dcf'
      INCLUDE 'incl/aerosol_component.dcf'
      INCLUDE 'incl/input_head.dcf'
!
!
!
!     SET THE COMPONENTS OF THE STANDARD PROFILES.
      DATA
     &     TYPE_AEROSOL(1)/IP_WATER_SOLUBLE/
     &   , TYPE_AEROSOL(2)/IP_DUST_LIKE/
     &   , TYPE_AEROSOL(3)/IP_OCEANIC/
     &   , TYPE_AEROSOL(4)/IP_SOOT/
     &   , TYPE_AEROSOL(5)/IP_ASH/
     &   , TYPE_AEROSOL(6)/IP_SULPHURIC/
!
!
!     OBTAIN THE FILES CONTAINING THE GRID ON WHICH TO SET THE PROFILE.
      WRITE(IU_USER, '(/A)')
     &   'ENTER THE BASE-NAME OF THE FILES TO BE USED.'
      READ(IU_STDIN, '(A)') BASE_NAME
      J=LEN(BASE_NAME)
1     IF (BASE_NAME(J:J).EQ.' ') THEN
         J=J-1
         IF (J.EQ.0) THEN
            WRITE(IU_ERR, '(A)') 'NO NAME SUPPLIED.'
            IERR=I_ERR_FATAL
            STOP
         ELSE
            GOTO 1
         ENDIF
      ENDIF
      LENGTH_NAME=J
!
!     READ THE SURFACE PRESSURE, THE TEMPERATURES AND SPECIFIC HUMIDITY.
      FILE_NAME(1: LENGTH_NAME+7)=BASE_NAME(1: LENGTH_NAME)
     &   //'.'//PHYS_SUFFIX(IP_PRESSURE_GROUND)
      CALL READ_SAMSON_FIELD(IERR, .TRUE., .FALSE.
     &   , FILE_NAME, TITLE_DUMMY
     &   , NDUM, N_LAYER_DUMMY, P_DUMMY, PSTAR_DUMMY, PSTAR)
      IF (IERR.EQ.I_MISSING_DATA) THEN
         WRITE(IU_ERR, '(/, 3A)') '*** ERROR: FILE ', FILE_NAME
     &      , ' DOES NOT EXIST'
         IERR=I_ERR_FATAL
         STOP
      ENDIF
      FILE_NAME(1: LENGTH_NAME+7)=BASE_NAME(1: LENGTH_NAME)
     &   //'.'//PHYS_SUFFIX(IP_TEMPERATURE)
      CALL READ_SAMSON_FIELD(IERR, .TRUE., .FALSE.
     &   , FILE_NAME, TITLE_DUMMY
     &   , NDUM, N_LAYER, P, PSTAR, T)
      IF (IERR.EQ.I_MISSING_DATA) THEN
         WRITE(IU_ERR, '(/, 3A)') '*** ERROR: FILE ', FILE_NAME
     &      , ' DOES NOT EXIST'
         IERR=I_ERR_FATAL
         STOP
      ENDIF
!     READ THE SPECIFIC HUMIDITY
      FILE_NAME(1: LENGTH_NAME+7)=BASE_NAME(1: LENGTH_NAME)
     &   //'.'//PHYS_SUFFIX(IP_SPEC_HUMIDITY)
      CALL READ_SAMSON_FIELD(IERR, .TRUE., .FALSE.
     &   , FILE_NAME, TITLE_DUMMY
     &   , NDUM, N_LAYER, P, PSTAR, Q)
      IF (IERR.EQ.I_MISSING_DATA) THEN
         WRITE(IU_ERR, '(/, 3A)') '*** ERROR: FILE ', FILE_NAME
     &      , ' DOES NOT EXIST'
         IERR=I_ERR_FATAL
         STOP
      ENDIF
!
!     INTEGRATE THROUGH THE ATMOSPHERE TO FIND THE HEIGHTS AT EACH
!     SPECIFIED LEVEL. THE SURFACE TEMPERATURE AND SPECIFIC HUMIDITY
!     ARE EXTRAPOLATED.
      TSTAR=(PSTAR*(T(N_LAYER)-T(N_LAYER-1))+P(N_LAYER)
     &   *T(N_LAYER-1)-P(N_LAYER-1)*T(N_LAYER))
     &   /(P(N_LAYER)-P(N_LAYER-1))
      QSTAR=(PSTAR*(Q(N_LAYER)-Q(N_LAYER-1))+P(N_LAYER)
     &   *Q(N_LAYER-1)-P(N_LAYER-1)*Q(N_LAYER))
     &   /(P(N_LAYER)-P(N_LAYER-1))
      DENSITY_STAR=FNC_DENSITY(PSTAR, TSTAR, QSTAR)
      DENSITY(N_LAYER)=FNC_DENSITY(P(N_LAYER), T(N_LAYER)
     &   , Q(N_LAYER))
      Z(N_LAYER)=(PSTAR-P(N_LAYER))/(GRAV_ACC*0.5D+00
     &   *(DENSITY_STAR+DENSITY(N_LAYER)))
      DO I=N_LAYER-1, 0, -1
         DENSITY(I)=FNC_DENSITY(P(I), T(I), Q(I))
         Z(I)=Z(I+1)+(P(I+1)-P(I))
     &      /(GRAV_ACC*0.5D+00*(DENSITY(I+1)+DENSITY(I)))
      ENDDO
!
!     ASSIGN THE CONCENTRATION VALUES FOR EACH TYPE OF AEROSOL.
      WRITE(IU_USER, '(/A)')
     &   'ENTER THE TYPE OF THE SRA-PROFILE TO BE SET.'
2     READ(IU_STDIN, *, IOSTAT=IOS) I_AEROSOL_PROFILE
      IF (IOS.NE.0) THEN
         WRITE(IU_ERR, '(A)') '+++ ILLEGAL RESPONSE:'
         IF (LOCK_CODE(.TRUE.)) THEN
            INCLUDE 'batch_error_main.finc'
         ELSE
            WRITE(IU_USER, '(A)') 'PLEASE RETYPE.'
            GOTO 2
         ENDIF
      ENDIF

      IF ( (I_AEROSOL_PROFILE.LE.0)
     &   .OR.(I_AEROSOL_PROFILE.GT.NPD_AEROSOL_PROFILE) ) THEN
         WRITE(IU_ERR, '(/A)') 
     &      '*** ERROR: UNRECOGINIZED AEROSOL PROFILE.'
         IF (LOCK_CODE(.TRUE.)) THEN
            INCLUDE 'batch_error_main.finc'
         ELSE
            WRITE(IU_USER, '(A)') 'PLEASE REENTER.'
            GOTO 2
         ENDIF
      ELSE
!
!        SET ALL VOLUME FRACTIONS TO ZERO. DO-LOOP PREFERRED
!        TO DATA STATEMENT IN CASE OF MULTIPLE CALLS.
         DO I=0, N_LAYER
            DO J=1, NPD_AEROSOL_COMPONENT
               VOLUME_FRACTION(I, J)=0.0D+00
            ENDDO
         ENDDO
!
         DO I=0, N_LAYER
!           FIND THE LAYER OF AEROSOL AT THIS HEIGHT
            I_AEROSOL_LAYER=1
3           IF (Z(I).GT.Z_AEROSOL_LAYER(I_AEROSOL_LAYER
     &         , I_AEROSOL_PROFILE)-TOL_TEST) THEN
               I_AEROSOL_LAYER=I_AEROSOL_LAYER+1
               IF (I_AEROSOL_LAYER.LT.
     &            N_AEROSOL_LAYER(I_AEROSOL_PROFILE)) GOTO 3
            ENDIF
!
            IF (I_AEROSOL_LAYER.LT.
     &         N_AEROSOL_LAYER(I_AEROSOL_PROFILE)) THEN
!              CALCULATE THE OVERALL EXTINCTION AT 550NM FOR THIS PROFILE.
!              ABOVE THE TOP OF THE PROFILES THE VOLUME FRACTION IS
!              LEFT AT 0.0.
               EXT_550NM_OVERALL=EXTINCTION_PROFILE(Z(I)
     &            , I_EXTINCTION_REP(I_AEROSOL_LAYER, I_AEROSOL_PROFILE)
     &            , NPD_EXTINCTION_PARAMETER
     &            , EXTINCTION_PARAMETER(1, I_AEROSOL_LAYER
     &            , I_AEROSOL_PROFILE)
     &            )
!
!              WORK OUT THE APPROPRIATE VOLUME FRACTIONS FOR THE MODEL
!              APPLICABLE IN THIS LAYER.
               I_MODEL=I_AEROSOL_MODEL(I_AEROSOL_LAYER
     &            , I_AEROSOL_PROFILE)
               CALL CALC_VOLUME_FRACTION(NPD_AEROSOL_COMPONENT
     &            , FRACTION_COMPONENT(1, I_MODEL)
     &            , EXT_550NM_COMPONENT, EXT_550NM_OVERALL
     &            , VOLUME_FRACTION_LEVEL
     &            )
               DO J=1, NPD_AEROSOL_COMPONENT
                  VOLUME_FRACTION(I, J)=VOLUME_FRACTION_LEVEL(J)
               ENDDO
            ENDIF
         ENDDO
!
!        CONVERT VOLUME FRACTION TO A MIXING RATIO.
         DO I=0, N_LAYER
            DO J=1, 6
               I_COMPONENT=TYPE_AEROSOL(J)
               AEROSOL_MIX_RATIO(I, J)
     &            =VOLUME_FRACTION(I, I_COMPONENT)
     &            *DENSITY_COMPONENT(I_COMPONENT)
     &            /DENSITY(I)
            ENDDO
         ENDDO
!
      ENDIF
!
!     WRITE OUT THE FILES OF AEROSOL COMPONENTS.
      DO J=1, 6
         I_COMPONENT=TYPE_AEROSOL(J)
         FILE_NAME(1: LENGTH_NAME+7)=BASE_NAME(1: LENGTH_NAME)
     &      //'.'//AEROSOL_SUFFIX(I_COMPONENT)
         CALL WRITE_SAMSON_FIELD(IERR, .FALSE.
     &      , FILE_NAME, AEROSOL_TITLE(I_COMPONENT)
     &      , 1, N_LAYER, P, PSTAR, AEROSOL_MIX_RATIO(0, I_COMPONENT))
      ENDDO

!
!
      STOP
      END
