!+ Program to run the radiation code with given optical properties. with 
!
! Method:
!	Files of optical properties are read in, options are selected
!	and the appropriate solver is called.
!
! 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 DRIVER_SOLVER
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_field.ucf'
      INCLUDE 'incl/dimensions_fixed.pcf'
      INCLUDE 'incl/dimensions_spec.ucf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/spectral_region.pcf'
      INCLUDE 'incl/cloud_scheme.pcf'
      INCLUDE 'incl/angular_integration.pcf'
      INCLUDE 'incl/solver.pcf'
      INCLUDE 'incl/gas_list.pcf'
      INCLUDE 'incl/aerosol_component.pcf'
      INCLUDE 'incl/input_head.pcf'
      INCLUDE 'incl/input_head.ncf'
      INCLUDE 'incl/physical_constants_0.ccf'
      INCLUDE 'incl/math_cnst.ccf'
      INCLUDE 'incl/two_stream_scheme.pcf'
      INCLUDE 'incl/prec_machine.ucf'
      INCLUDE 'incl/precision.pcf'
      INCLUDE 'incl/error.pcf'
!
!     DECLARATION OF VARIABLES.
      CHARACTER
     &     BASE_NAME_FREE*80
     &   , BASE_NAME_CLOUD*80
     &   , FILE_NAME*80
     &   , TITLE_DUMMY*80
     &   , CHAR_YN*1
!
      LOGICAL
     &     LOCK_CODE
!             LOGICAL TO FORBID INTERACTIVE LOOPING
      EXTERNAL
     &     LOCK_CODE
!
      LOGICAL
     &     L_LEVEL_SET
     &   , L_CLOUD
     &   , L_2_STREAM_CORRECT
     &   , L_RESCALE
     &   , L_NET
     &   , L_DUMMY
      INTEGER
     &     IERR
     &   , N_LAYER
     &   , N_LAYER_DUMMY
     &   , N_AUGMENT
     &   , ISOLIR
     &   , N_FREE_PROFILE(NPD_LAYER)
     &   , I_FREE_PROFILE(NPD_PROFILE, NPD_LAYER)
     &   , N_CLOUD_PROFILE(NPD_LAYER)
     &   , I_CLOUD_PROFILE(NPD_PROFILE, NPD_LAYER)
     &   , I_CLOUD
     &   , I_ANGULAR_INTEGRATION
     &   , I_SOLVER
!             MAIN SOLVER
     &   , I_EXTRAPOLATE_MIX_NET
!             EXTRAPOLATION SCHEME FOR THE MIXED SOLVER
     &   , I_2STREAM
     &   , N_ORDER_GAUSS
     &   , J
     &   , I
     &   , LENGTH_NAME_FREE
     &   , LENGTH_NAME_CLOUD
     &   , I_DUMMY
      INTEGER
     &     NDUM
      REAL*8
     &     TAU_FREE(NPD_PROFILE, NPD_LAYER)
     &   , OMEGA_FREE(NPD_PROFILE, NPD_LAYER)
     &   , ASYMMETRY_FREE(NPD_PROFILE, NPD_LAYER)
     &   , FORWARD_SCATTER_FREE(NPD_PROFILE, NPD_LAYER)
      REAL*8
     &     TAU_CLOUD(NPD_PROFILE, NPD_LAYER)
     &   , OMEGA_CLOUD(NPD_PROFILE, NPD_LAYER)
     &   , ASYMMETRY_CLOUD(NPD_PROFILE, NPD_LAYER)
     &   , FORWARD_SCATTER_CLOUD(NPD_PROFILE, NPD_LAYER)
      REAL*8
     &     SEC_0(NPD_PROFILE)
     &   , ALBEDO_SURFACE_DIR(NPD_PROFILE)
     &   , ALBEDO_SURFACE_DIFF(NPD_PROFILE)
     &   , GROUND_EMISSION(NPD_PROFILE)
!             FLUX EMITTED FROM GROUND
     &   , SOURCE_GROUND(NPD_PROFILE)
!             SOURCE TERM AT GROUND
     &   , FLUX_INC_DOWN(NPD_PROFILE)
     &   , FLUX_INC_DIRECT(NPD_PROFILE)
      REAL*8
     &     P(NPD_PROFILE, 0: NPD_LAYER)
     &   , W_CLOUD(NPD_PROFILE, NPD_LAYER)
     &   , FRAC_CLOUD(NPD_PROFILE, NPD_LAYER, NPD_CLOUD_TYPE)
!             FRACTIONS OF TYPES OF CLOUDS
     &   , PLANCK_SOURCE(NPD_PROFILE, 0: NPD_LAYER)
     &   , DIFF_PLANCK(NPD_PROFILE, NPD_LAYER)
     &   , R_DUMMY(NPD_PROFILE, NPD_LAYER)
      REAL*8
     &     FLUX_DIRECT(NPD_PROFILE, 0: NPD_LAYER)
!             DIRECT FLUX
     &   , FLUX_DIFFUSE(NPD_PROFILE, 2*NPD_LAYER+2)
!             DIFFUSE FLUX
     &   , FLUX_DOWN(NPD_PROFILE, 0: NPD_LAYER)
!             DOWNWARD FLUX
     &   , FLUX_TOTAL_DOWN(NPD_PROFILE, 0: NPD_LAYER)
!             TOTAL DOWNWARD FLUX
     &   , FLUX_UP(NPD_PROFILE, 0: NPD_LAYER)
!             UPWARD FLUX
     &   , FLUX_NET(NPD_PROFILE, 0: NPD_LAYER)
!             NET FLUX
     &   , FLUX_DIRECT_DUMMY(NPD_PROFILE, 0: NPD_LAYER)
!             DIRECT FLUX
     &   , FLUX_DIFFUSE_DUMMY(NPD_PROFILE, 2*NPD_LAYER+2)
!             DIFFUSE FLUX
      REAL*8
     &     PSTAR(NPD_PROFILE)
      REAL*8
     &     W_FREE(NPD_PROFILE, NPD_LAYER)
!             CLOUD-FREE AMOUNTS
     &   , CLOUD_OVERLAP(NPD_PROFILE, 0: NPD_LAYER, NPD_OVERLAP_COEFF)
!             COEFFICIENTS FOR TRANSFER OF ENERGY AT INTERFACE
      INTEGER
     &     N_COLUMN(NPD_PROFILE)
!             NUMBER OF COLUMNS REQUIRED
      LOGICAL
     &     L_COLUMN(NPD_PROFILE, NPD_LAYER, NPD_COLUMN)
!             COLUMN FLAGS FOR COLUMNS
      REAL*8
     &     AREA_COLUMN(NPD_PROFILE, NPD_COLUMN)
!             AREAS OF COLUMNS
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     READ_SAMSON_FIELD, OVERLAP_MIX_MAXIMUM
     &   , OVERLAP_MIX_RANDOM, SPLIT_MAXIMUM, RESCALE_ASYMMETRY
     &   , RESCALE_TAU_OMEGA
     &   , TWO_STREAM, MIX_COLUMN, CLOUD_COLUMN
     &   , GAUSS_ANGLE, WRITE_SAMSON_FIELD
!
!
      INCLUDE 'incl/input_head.dcf'
      INCLUDE 'incl/samson_pstar.ecf'
!
      DATA IERR/I_NORMAL/
      DATA PSTAR/NPD_PROFILE*PSTAR_DEFAULT/
!
!
!
      WRITE(IU_USER, '(/A)') 'ENTER TYPE NUMBER OF CLOUD SCHEME.'
      READ(IU_STDIN, *) I_CLOUD
!     SET THE CLOUD FLAG FROM CONSIDERATION OF THE VALUE OF I_CLOUD.
      IF ( (I_CLOUD.EQ.IP_CLOUD_MIX_MAX).OR.
     &     (I_CLOUD.EQ.IP_CLOUD_MIX_RANDOM).OR.
     &     (I_CLOUD.EQ.IP_CLOUD_COLUMN_MAX) ) THEN
!        CLOUDS ARE INCLUDED.
         L_CLOUD=.TRUE.
      ELSE
!        CLOUDS ARE NOT INCLUDED.
         L_CLOUD=.FALSE.
      ENDIF
!
      WRITE(IU_USER, '(/A)') 'ENTER TYPE OF ANGULAR INTEGRATION.'
      READ(IU_STDIN, *) I_ANGULAR_INTEGRATION
      IF (I_CLOUD.NE.IP_IR_GAUSS) THEN
         WRITE(IU_USER, '(A)') 'ENTER TWO-STREAM SCHEME.'
         READ(IU_STDIN, *) I_2STREAM
         WRITE(IU_USER, '(/A)') 'ENTER RESCALING FLAG. (T/F)'
         READ(IU_STDIN, *) L_RESCALE
         IF (ISOLIR.EQ.IP_INFRA_RED) THEN
            WRITE (IU_USER, '(/A)')
     &         'DO YOU WANT TO APPLY THE EDGE CORRECTION? (Y/N)'
            READ(IU_STDIN, '(A)') CHAR_YN
            IF ( (CHAR_YN.EQ.'Y').OR.(CHAR_YN.EQ.'y') ) THEN
               L_2_STREAM_CORRECT=.TRUE.
            ELSE
               L_2_STREAM_CORRECT=.FALSE.
            ENDIF
         ENDIF
         WRITE(IU_USER, '(/A)') 'ENTER TYPE NUMBER OF SOLVER.'
         READ(IU_STDIN, *) I_SOLVER
!
!        SET THE SOLUTION FOR NET FLUXES.
         IF ( (I_SOLVER.EQ.IP_SOLVER_TRIDIAG_NET).OR.
     &        (I_SOLVER.EQ.IP_SOLVER_MIX_NET) ) THEN
            L_NET=.TRUE.
         ELSE
            L_NET=.FALSE.
         ENDIF
!
         IF (I_SOLVER.EQ.IP_SOLVER_MIX_NET) THEN
            WRITE(IU_USER, '(A)') 'EXTRAPOLATION SCHEME FOR SOLVER.'
            READ(IU_STDIN, *) I_EXTRAPOLATE_MIX_NET
         ENDIF
!
      ELSE IF (I_CLOUD.EQ.IP_IR_GAUSS) THEN
!
         WRITE(IU_USER, '(/A)')
     &      'ENTER THE ORDER OF GAUSSIAN QUADRATURE.'
         READ(IU_STDIN, *) N_ORDER_GAUSS
         L_NET=.FALSE.
!
      ENDIF
!
!
!
!     READ INPUT FILES.
      WRITE(IU_USER, '(A)') 'ENTER BASE NAME OF CLEAR-SKY INPUT FILES.'
      READ(IU_STDIN, '(A)') BASE_NAME_FREE
      J=LEN(BASE_NAME_FREE)
1     IF (BASE_NAME_FREE(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_FREE=J
!
!
!     THE FIELDS HERE ARE GIVEN IN LAYERS, BUT HAVE NO 0 INDEXING AT THE
!     TOP: THEY MUST THEREFORE BE READ WITHOUT AN OFFSET AND THE NUMBER
!     OF LAYERS MUST BE ADJUSTED.
!
      FILE_NAME(1: LENGTH_NAME_FREE+7)
     &   =BASE_NAME_FREE(1: LENGTH_NAME_FREE)
     &   //'.'//PHYS_SUFFIX(IP_TAU)
      CALL READ_SAMSON_FIELD(IERR, .TRUE., .FALSE.
     &   , FILE_NAME, TITLE_DUMMY
     &   , NDUM, N_LAYER_DUMMY, P, PSTAR, TAU_FREE)
      N_LAYER=N_LAYER_DUMMY+1
      IF (IERR.EQ.I_MISSING_DATA) THEN
         WRITE(IU_ERR, '(/A)')
     &      '*** ERROR IN PROGRAM SAMSON_DRIVER.'
         WRITE(IU_ERR, '(3A)') 'FILE ', FILE_NAME
     &      , ' DOES NOT EXIST'
         IERR=I_ERR_FATAL
         STOP
      ELSE
         L_LEVEL_SET=.TRUE.
      ENDIF
!
      FILE_NAME(1: LENGTH_NAME_FREE+7)
     &   =BASE_NAME_FREE(1: LENGTH_NAME_FREE)
     &   //'.'//PHYS_SUFFIX(IP_SSA)
      CALL READ_SAMSON_FIELD(IERR, .TRUE., .FALSE.
     &   , FILE_NAME, TITLE_DUMMY
     &   , NDUM, N_LAYER_DUMMY, P, PSTAR, OMEGA_FREE)
      IF (IERR.EQ.I_MISSING_DATA) THEN
         WRITE(IU_ERR, '(/A)')
     &      '*** ERROR IN PROGRAM SAMSON_DRIVER.'
         WRITE(IU_ERR, '(3A)') 'FILE ', FILE_NAME
     &      , ' DOES NOT EXIST'
         IERR=I_ERR_FATAL
         STOP
      ELSE
         L_LEVEL_SET=.TRUE.
      ENDIF
!
      FILE_NAME(1: LENGTH_NAME_FREE+7)
     &   =BASE_NAME_FREE(1: LENGTH_NAME_FREE)
     &   //'.'//PHYS_SUFFIX(IP_GSC)
      CALL READ_SAMSON_FIELD(IERR, .TRUE., .FALSE.
     &   , FILE_NAME, TITLE_DUMMY
     &   , NDUM, N_LAYER_DUMMY, P, PSTAR, ASYMMETRY_FREE)
      IF (IERR.EQ.I_MISSING_DATA) THEN
         WRITE(IU_ERR, '(/A)')
     &      '*** ERROR IN PROGRAM SAMSON_DRIVER.'
         WRITE(IU_ERR, '(3A)') 'FILE ', FILE_NAME
     &      , ' DOES NOT EXIST'
         IERR=I_ERR_FATAL
         STOP
      ELSE
         L_LEVEL_SET=.TRUE.
      ENDIF
!
      IF (L_RESCALE) THEN
         FILE_NAME(1: LENGTH_NAME_FREE+7)
     &      =BASE_NAME_FREE(1: LENGTH_NAME_FREE)
     &      //'.'//PHYS_SUFFIX(IP_FSC)
         CALL READ_SAMSON_FIELD(IERR, .TRUE., .FALSE.
     &      , FILE_NAME, TITLE_DUMMY
     &      , NDUM, N_LAYER_DUMMY, P, PSTAR, FORWARD_SCATTER_FREE)
         IF (IERR.EQ.I_MISSING_DATA) THEN
            WRITE(IU_ERR, '(/A)')
     &         '*** ERROR IN PROGRAM SAMSON_DRIVER.'
            WRITE(IU_ERR, '(3A)') 'FILE ', FILE_NAME
     &         , ' DOES NOT EXIST'
            IERR=I_ERR_FATAL
            STOP
         ELSE
            L_LEVEL_SET=.TRUE.
         ENDIF
      ENDIF
!
!
!     READ IN THE CLOUD DATA WHEN CLOUDS ARE PRESENT.
      IF (L_CLOUD) THEN
         WRITE(IU_USER, '(A)') 'ENTER BASE NAME OF CLOUDY INPUT FILES.'
         READ(IU_STDIN, '(A)') BASE_NAME_CLOUD
         J=LEN(BASE_NAME_CLOUD)
2        IF (BASE_NAME_CLOUD(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 2
            ENDIF
         ENDIF
         LENGTH_NAME_CLOUD=J
!
         FILE_NAME(1: LENGTH_NAME_CLOUD+7)
     &      =BASE_NAME_CLOUD(1: LENGTH_NAME_CLOUD)
     &      //'.'//PHYS_SUFFIX(IP_TAU)
         CALL READ_SAMSON_FIELD(IERR, .TRUE., .FALSE.
     &      , FILE_NAME, TITLE_DUMMY
     &      , NDUM, N_LAYER_DUMMY, P, PSTAR, TAU_CLOUD)
         IF (IERR.EQ.I_MISSING_DATA) THEN
            WRITE(IU_ERR, '(/A)')
     &         '*** ERROR IN PROGRAM SAMSON_DRIVER.'
            WRITE(IU_ERR, '(3A)') 'FILE ', FILE_NAME
     &         , ' DOES NOT EXIST'
            IERR=I_ERR_FATAL
            STOP
         ELSE
            L_LEVEL_SET=.TRUE.
         ENDIF
!
         FILE_NAME(1: LENGTH_NAME_CLOUD+7)
     &      =BASE_NAME_CLOUD(1: LENGTH_NAME_CLOUD)
     &      //'.'//PHYS_SUFFIX(IP_SSA)
         CALL READ_SAMSON_FIELD(IERR, .TRUE., .FALSE.
     &      , FILE_NAME, TITLE_DUMMY
     &      , NDUM, N_LAYER_DUMMY, P, PSTAR, OMEGA_CLOUD)
         IF (IERR.EQ.I_MISSING_DATA) THEN
            WRITE(IU_ERR, '(/A)')
     &         '*** ERROR IN PROGRAM SAMSON_DRIVER.'
            WRITE(IU_ERR, '(3A)') 'FILE ', FILE_NAME
     &         , ' DOES NOT EXIST'
            IERR=I_ERR_FATAL
            STOP
         ELSE
            L_LEVEL_SET=.TRUE.
         ENDIF
!
         FILE_NAME(1: LENGTH_NAME_CLOUD+7)
     &      =BASE_NAME_CLOUD(1: LENGTH_NAME_CLOUD)
     &      //'.'//PHYS_SUFFIX(IP_GSC)
         CALL READ_SAMSON_FIELD(IERR, .TRUE., .FALSE.
     &      , FILE_NAME, TITLE_DUMMY
     &      , NDUM, N_LAYER_DUMMY, P, PSTAR, ASYMMETRY_CLOUD)
         IF (IERR.EQ.I_MISSING_DATA) THEN
            WRITE(IU_ERR, '(/A)')
     &         '*** ERROR IN PROGRAM SAMSON_DRIVER.'
            WRITE(IU_ERR, '(3A)') 'FILE ', FILE_NAME
     &         , ' DOES NOT EXIST'
            IERR=I_ERR_FATAL
            STOP
         ELSE
            L_LEVEL_SET=.TRUE.
         ENDIF
!
         IF (L_RESCALE) THEN
            FILE_NAME(1: LENGTH_NAME_CLOUD+7)
     &         =BASE_NAME_CLOUD(1: LENGTH_NAME_CLOUD)
     &         //'.'//PHYS_SUFFIX(IP_FSC)
            CALL READ_SAMSON_FIELD(IERR, .TRUE., .FALSE.
     &         , FILE_NAME, TITLE_DUMMY
     &         , NDUM, N_LAYER_DUMMY, P, PSTAR, FORWARD_SCATTER_CLOUD)
            IF (IERR.EQ.I_MISSING_DATA) THEN
               WRITE(IU_ERR, '(/A)')
     &            '*** ERROR IN PROGRAM SAMSON_DRIVER.'
               WRITE(IU_ERR, '(3A)') 'FILE ', FILE_NAME
     &            , ' DOES NOT EXIST'
               IERR=I_ERR_FATAL
               STOP
            ELSE
               L_LEVEL_SET=.TRUE.
            ENDIF
         ENDIF
!
         FILE_NAME(1: LENGTH_NAME_CLOUD+7)
     &      =BASE_NAME_CLOUD(1: LENGTH_NAME_CLOUD)
     &      //'.'//PHYS_SUFFIX(IP_CLOUD_FRACTION)
         CALL READ_SAMSON_FIELD(IERR, .TRUE., .FALSE.
     &      , FILE_NAME, TITLE_DUMMY
     &      , NDUM, N_LAYER_DUMMY, P, PSTAR, W_CLOUD)
         IF (IERR.EQ.I_MISSING_DATA) THEN
            WRITE(IU_ERR, '(/A)')
     &         '*** ERROR IN PROGRAM SAMSON_DRIVER.'
            WRITE(IU_ERR, '(3A)') 'FILE ', FILE_NAME
     &         , ' DOES NOT EXIST'
            IERR=I_ERR_FATAL
            STOP
         ELSE
            L_LEVEL_SET=.TRUE.
         ENDIF
!
!        SET THE CLOUDS IN THE LAYERS.
         DO I=1, N_LAYER
            N_FREE_PROFILE(I)=0
            N_CLOUD_PROFILE(I)=0
            FRAC_CLOUD(1, I, 1)=1.0D+00
            IF (W_CLOUD(1, I).GT.TOL_TEST) THEN
               N_CLOUD_PROFILE(I)=1
               I_CLOUD_PROFILE(1, I)=1
            ELSE
               N_FREE_PROFILE(I)=1
               I_FREE_PROFILE(1, I)=1
            ENDIF
         ENDDO
!
      ENDIF
!
!     OBTAIN THE BOUNDARY CONDITIONS.
      WRITE(IU_USER, '(/A)') 'ENTER SPECTRAL REGION.'
3     READ(IU_STDIN, *) ISOLIR
      IF (ISOLIR.EQ.IP_SOLAR) THEN
         WRITE(IU_USER, '(A)') 'ENTER SOLAR ZENITH ANGLE IN DEGREES.'
         READ(IU_STDIN, *) SEC_0(1)
         SEC_0(1)=1.0D+00/COS(SEC_0(1)*PI/180.0D+00)
         WRITE(IU_USER, '(A)') 'ENTER SOLAR IRRADIANCE AT TOA.'
         READ(IU_STDIN, *) FLUX_INC_DIRECT(1)
         FLUX_INC_DIRECT(1)=FLUX_INC_DIRECT(1)/SEC_0(1)
         WRITE(IU_USER, '(A)') 'ENTER SURFACE ALBEDO FOR SOLAR BEAM.'
         READ(IU_STDIN, *) ALBEDO_SURFACE_DIR(1)
      ELSE IF (ISOLIR.EQ.IP_INFRA_RED) THEN
         WRITE(IU_USER, '(A)') 'ENTER SOURCE FUNCTION OF GROUND.'
         READ(IU_STDIN, *) GROUND_EMISSION(1)
      ELSE
         WRITE(IU_ERR, '(A)') '+++ UNKNOWN SPECTRAL REGION:'
         IF (LOCK_CODE(.TRUE.)) THEN
            INCLUDE 'batch_error_main.finc'
         ELSE
            WRITE(IU_USER, '(A)') 'PLEASE REENTER.'
            GOTO 3
         ENDIF
      ENDIF
      WRITE(IU_USER, '(A)') 'ENTER SURFACE ALBEDO FOR DIFFUSE BEAMS.'
      READ(IU_STDIN, *) ALBEDO_SURFACE_DIFF(1)
!     RESET THE DIFFUSE ALBEDO TO AVOID NUMERICAL INSTABILITY.
      ALBEDO_SURFACE_DIFF(1)=MAX(TOL_DIV, ALBEDO_SURFACE_DIFF(1))
      WRITE(IU_USER, '(A)') 'ENTER INCIDENT DIFFUSE FLUX.'
      READ(IU_STDIN, *) FLUX_INC_DOWN(1)
!
!     IN THE INFRA_RED THE PLANCK FUNCTION IS ALSO NEEDED.
      IF (ISOLIR.EQ.IP_INFRA_RED) THEN
         FILE_NAME(1: LENGTH_NAME_FREE+7)
     &      =BASE_NAME_FREE(1: LENGTH_NAME_FREE)
     &      //'.'//PHYS_SUFFIX(IP_PLANCK)
         CALL READ_SAMSON_FIELD(IERR, .TRUE., .TRUE.
     &      , FILE_NAME, TITLE_DUMMY
     &      , NDUM, N_LAYER_DUMMY, P, PSTAR, PLANCK_SOURCE)
         IF (IERR.EQ.I_MISSING_DATA) THEN
            WRITE(IU_ERR, '(/A)')
     &         '*** ERROR IN PROGRAM SAMSON_DRIVER.'
            WRITE(IU_ERR, '(3A)') 'FILE ', FILE_NAME
     &         , ' DOES NOT EXIST'
            IERR=I_ERR_FATAL
            STOP
         ELSE
            L_LEVEL_SET=.TRUE.
         ENDIF
      ELSE IF (ISOLIR.EQ.IP_SOLAR) THEN
!        CONVERT TO TOTAL DOWNWARD FLUX.
         FLUX_INC_DOWN(1)=FLUX_INC_DOWN(1)+FLUX_INC_DIRECT(1)
      ENDIF
!
!     INPUT OF DATA COMPLETED. PROCEED TO SOLVE EQUATIONS.
!
!     THE LENGTH OF THE LONG FLUX VECTOR DEPENDS ON WHETHER WE SOLVE
!     FOR THE FULL OR THE NET FLUX.
      IF (L_NET) THEN
         N_AUGMENT=N_LAYER+1
      ELSE
         N_AUGMENT=2*(N_LAYER+1)
      ENDIF
!
!     SET THE GEOMETRY FOR THE CLOUD ALGORITHM
      IF (I_CLOUD.EQ.IP_CLOUD_MIX_MAX) THEN
         CALL OVERLAP_MIX_MAXIMUM(1, N_LAYER, 1
     &      , ISOLIR, I_SOLVER
     &      , W_CLOUD, W_FREE
     &      , CLOUD_OVERLAP
     &      )
      ELSE IF (I_CLOUD.EQ.IP_CLOUD_MIX_RANDOM) THEN
         CALL OVERLAP_MIX_RANDOM(1, N_LAYER, 1
     &      , ISOLIR, I_SOLVER
     &      , W_CLOUD, W_FREE
     &      , CLOUD_OVERLAP
     &      )
      ELSE IF (I_CLOUD.EQ.IP_CLOUD_COLUMN_MAX) THEN
!        SPLIT THE ATMOSPHERE INTO COLUMNS.
         CALL SPLIT_MAXIMUM(1, N_LAYER
     &      , W_CLOUD
     &      , N_COLUMN, AREA_COLUMN, L_COLUMN
     &      )
      ENDIF
!
!     THE DIFFERENCES IN THE THERMAL SOURCE FUNCTION ARE CALCULATED
!     DIRECTLY HERE, AND NOT BY THE STANDARD ROUTINE SINCE THEY ARE
!     ARE NOT HERE BASED ON A POLYNOMIAL FIT IN THE TEMPERATURE.
      IF (ISOLIR.EQ.IP_INFRA_RED) THEN
         DO I=1, N_LAYER
            DIFF_PLANCK(1, I)
     &         =PLANCK_SOURCE(1, I)-PLANCK_SOURCE(1, I-1)
         ENDDO
!        SINCE THE FLUXES CALCULATED ARE DIFFERENTIAL THE SOURCES
!        AT THE BOUNDARIES MUST BE ALTERED.
         FLUX_INC_DOWN(1)=FLUX_INC_DOWN(1)-PLANCK_SOURCE(1, 0)
         SOURCE_GROUND(1)=GROUND_EMISSION(1)-PLANCK_SOURCE(1, N_LAYER)
      ENDIF
!
!     FULL CLOUDY CALCULATIONS ARE AVAILABLE ONLY WITH THE
!     TWO-STREAM APPROXIMATION AT PRESENT.
      IF (I_ANGULAR_INTEGRATION.EQ.IP_TWO_STREAM) THEN
!        CALCULATE THE SCATTERING COEFFICIENTS FROM THE INPUT DATA.
         IF (L_RESCALE) THEN
            CALL RESCALE_ASYMMETRY(1, 1, N_LAYER
     &         , ASYMMETRY_FREE, FORWARD_SCATTER_FREE
     &         )
            CALL RESCALE_TAU_OMEGA(1, 1, N_LAYER
     &         , TAU_FREE, OMEGA_FREE, FORWARD_SCATTER_FREE
     &         )
            IF (L_CLOUD) THEN
               CALL RESCALE_ASYMMETRY(1, 1, N_LAYER
     &            , ASYMMETRY_CLOUD, FORWARD_SCATTER_CLOUD
     &            )
               CALL RESCALE_TAU_OMEGA(1, 1, N_LAYER
     &            , TAU_CLOUD, OMEGA_CLOUD, FORWARD_SCATTER_CLOUD
     &            )
            ENDIF
         ENDIF
!
!
!        SOLVE THE EQUATIONS USING THE SCHEME INDICATED BY THE VALUE
!        OF I_CLOUD.
         IF (I_CLOUD.EQ.IP_CLOUD_CLEAR) THEN
!           A TWO-STREAM SCHEME WITH NO CLOUDS.
            CALL TWO_STREAM(IERR
     &         , 1, N_LAYER
     &         , I_2STREAM
     &         , L_2_STREAM_CORRECT, PLANCK_SOURCE, GROUND_EMISSION
     &         , L_NET, I_SOLVER
     &         , .FALSE., R_DUMMY
     &         , ISOLIR
     &         , DIFF_PLANCK
     &         , .FALSE., R_DUMMY
     &         , FLUX_INC_DOWN, FLUX_INC_DIRECT, SEC_0
     &         , ALBEDO_SURFACE_DIFF, ALBEDO_SURFACE_DIR
     &         , SOURCE_GROUND
     &         , TAU_FREE, OMEGA_FREE, ASYMMETRY_FREE
     &         , FLUX_DIRECT, FLUX_DIFFUSE
     &         , L_DUMMY, FLUX_DIRECT_DUMMY, FLUX_DIFFUSE_DUMMY
     &         )
!
         ELSEIF ( (I_CLOUD.EQ.IP_CLOUD_MIX_MAX)
     &          .OR. (I_CLOUD.EQ.IP_CLOUD_MIX_RANDOM) ) THEN
!           CLOUDS ARE TREATED USING ZDUNKOWSKI'S MIXED-COLUMN SCHEME.
!           THE GEOMETRY HAS BEEN SET BEFORE.
            CALL MIX_COLUMN(IERR
     &         , 1, N_LAYER
     &         , I_2STREAM
     &         , L_2_STREAM_CORRECT, PLANCK_SOURCE, GROUND_EMISSION
     &         , I_SOLVER, I_EXTRAPOLATE_MIX_NET, L_NET
     &         , .FALSE., R_DUMMY
     &         , ISOLIR
     &         , DIFF_PLANCK
     &         , .FALSE., R_DUMMY
     &         , FLUX_INC_DOWN, FLUX_INC_DIRECT, SEC_0
     &         , ALBEDO_SURFACE_DIFF, ALBEDO_SURFACE_DIR, SOURCE_GROUND
     &         , TAU_FREE, OMEGA_FREE, ASYMMETRY_FREE
     &         , 1, 1, FRAC_CLOUD
     &         , W_FREE, N_FREE_PROFILE, I_FREE_PROFILE
     &         , W_CLOUD, N_CLOUD_PROFILE, I_CLOUD_PROFILE
     &         , CLOUD_OVERLAP
     &         , TAU_CLOUD, OMEGA_CLOUD, ASYMMETRY_CLOUD
     &         , FLUX_DIRECT, FLUX_DIFFUSE
     &         , .FALSE., I_DUMMY
     &         , FLUX_DIRECT_DUMMY, FLUX_DIFFUSE_DUMMY
     &         )
         ELSEIF (I_CLOUD.EQ.IP_CLOUD_COLUMN_MAX) THEN
!           CLOUDS ARE TREATED ON THE ASSUMPTION OF MAXIMUM OVERLAP
!           IN A COLUMN MODEL.
            CALL CLOUD_COLUMN(IERR
     &         , 1, N_LAYER
     &         , I_2STREAM
     &         , L_2_STREAM_CORRECT, PLANCK_SOURCE, GROUND_EMISSION
     &         , I_SOLVER, L_NET, N_AUGMENT
     &         , .FALSE., R_DUMMY
     &         , ISOLIR
     &         , DIFF_PLANCK
     &         , .FALSE., R_DUMMY
     &         , FLUX_INC_DOWN, FLUX_INC_DIRECT, SEC_0
     &         , ALBEDO_SURFACE_DIFF, ALBEDO_SURFACE_DIR, SOURCE_GROUND
     &         , TAU_FREE, OMEGA_FREE, ASYMMETRY_FREE
     &         , 1, 1, FRAC_CLOUD
     &         , N_COLUMN, L_COLUMN, AREA_COLUMN
     &         , TAU_CLOUD, OMEGA_CLOUD, ASYMMETRY_CLOUD
     &         , FLUX_DIRECT, FLUX_DIFFUSE
     &         , .FALSE., I_DUMMY
     &         , FLUX_DIRECT_DUMMY, FLUX_DIFFUSE_DUMMY
     &         )
         ENDIF
!
      ELSE IF (I_ANGULAR_INTEGRATION.EQ.IP_IR_GAUSS) THEN
!        FULL ANGULAR RESOLUTION USING GASUSSIAN INTEGRATION.
         CALL GAUSS_ANGLE(IERR
     &      , 1, N_LAYER, L_NET, N_AUGMENT
     &      , N_ORDER_GAUSS
     &      , TAU_FREE
     &      , FLUX_INC_DOWN
     &      , DIFF_PLANCK, SOURCE_GROUND, ALBEDO_SURFACE_DIFF
     &      , FLUX_DIFFUSE
     &      , .FALSE., R_DUMMY
     &      )
      ENDIF
      IF (IERR.NE.I_NORMAL) STOP
!
!     WRITE OUT THE SAMSON FILES.
      IF (L_NET) THEN
         DO I=0, N_LAYER
            FLUX_NET(1, I)=FLUX_DIFFUSE(1, I+1)
         ENDDO
         FILE_NAME(1: LENGTH_NAME_FREE+7)
     &      =BASE_NAME_FREE(1: LENGTH_NAME_FREE)
     &      //'.'//PHYS_SUFFIX(IP_FLUX_NET)
         CALL WRITE_SAMSON_FIELD(IERR, .FALSE.
     &      , FILE_NAME, PHYS_TITLE(IP_FLUX_NET)
     &      , 1, N_LAYER, P, PSTAR, FLUX_NET)
      ELSE
         DO I=0, N_LAYER
            FLUX_UP(1, I)=FLUX_DIFFUSE(1, 2*I+1)
            FLUX_DOWN(1, I)=FLUX_DIFFUSE(1, 2*I+2)-FLUX_DIRECT(1, I)
            FLUX_TOTAL_DOWN(1, I)=FLUX_DIFFUSE(1, 2*I+2)
            FLUX_NET(1, I)=FLUX_DIRECT(1, I)+FLUX_DOWN(1, I)
     &         -FLUX_UP(1, I)
         ENDDO
         IF (ISOLIR.EQ.IP_INFRA_RED) THEN
!           MOVE FROM DIFFERENTIAL TO ACTUAL FLUXES.
            DO I=0, N_LAYER
               FLUX_UP(1, I)=FLUX_UP(1, I)+PLANCK_SOURCE(1, I)
               FLUX_DOWN(1, I)=FLUX_DOWN(1, I)+PLANCK_SOURCE(1, I)
            ENDDO
         ENDIF
         FILE_NAME(1: LENGTH_NAME_FREE+7)
     &      =BASE_NAME_FREE(1: LENGTH_NAME_FREE)
     &      //'.'//PHYS_SUFFIX(IP_FLUX_DOWN)
         CALL WRITE_SAMSON_FIELD(IERR, .FALSE.
     &      , FILE_NAME, PHYS_TITLE(IP_FLUX_DOWN)
     &      , 1, N_LAYER, P, PSTAR, FLUX_DOWN)
         FILE_NAME(1: LENGTH_NAME_FREE+7)
     &      =BASE_NAME_FREE(1: LENGTH_NAME_FREE)
     &      //'.'//PHYS_SUFFIX(IP_FLUX_UP)
         CALL WRITE_SAMSON_FIELD(IERR, .FALSE.
     &      , FILE_NAME, PHYS_TITLE(IP_FLUX_UP)
     &      , 1, N_LAYER, P, PSTAR, FLUX_UP)
         FILE_NAME(1: LENGTH_NAME_FREE+7)
     &      =BASE_NAME_FREE(1: LENGTH_NAME_FREE)
     &      //'.'//PHYS_SUFFIX(IP_FLUX_NET)
         CALL WRITE_SAMSON_FIELD(IERR, .FALSE.
     &      , FILE_NAME, PHYS_TITLE(IP_FLUX_NET)
     &      , 1, N_LAYER, P, PSTAR, FLUX_NET)
      ENDIF
      IF (ISOLIR.EQ.IP_SOLAR) THEN
         FILE_NAME(1: LENGTH_NAME_FREE+7)
     &      =BASE_NAME_FREE(1: LENGTH_NAME_FREE)
     &      //'.'//PHYS_SUFFIX(IP_FLUX_DIRECT)
         CALL WRITE_SAMSON_FIELD(IERR, .FALSE.
     &      , FILE_NAME, PHYS_TITLE(IP_FLUX_DIRECT)
     &      , 1, N_LAYER, P, PSTAR, FLUX_DIRECT)
         FILE_NAME(1: LENGTH_NAME_FREE+7)
     &      =BASE_NAME_FREE(1: LENGTH_NAME_FREE)
     &      //'.'//PHYS_SUFFIX(IP_FLUX_TOTAL_DOWN)
         CALL WRITE_SAMSON_FIELD(IERR, .FALSE.
     &      , FILE_NAME, PHYS_TITLE(IP_FLUX_TOTAL_DOWN)
     &      , 1, N_LAYER, P, PSTAR, FLUX_TOTAL_DOWN)
      ENDIF
!
!
      STOP
      END
