!+ Subroutine to treat the growth of particles with humidity.
!
!*********************
! ELLIE's VERSION WHICH DOES NOT ALLOW ALL THE AEROSOL COMPONENTS THAT 
! THE FULL VERSION 5 OF THIS DOES - I couldn't work out which needed to be
! changed 12/8/99
!
!
! Method:
!
!	THIS ROUTINE CHANGES THE LOGARITHMS OF THE MEAN AND STANDARD
!	DEVIATION OF A LOG-NORMAL DISTRIBUTION OF PARTICLE SIZES, TO
!	ACCOUNT FOR GROWTH DUE TO HUMIDITY.
!
!	AUTHOR: D.L.ROBERTS.
!
!	THIS IS VERSION 5 OF THE ROUTINE. (MAY 11TH 1995)
!	VERSION 1 DID NOT ALTER THE STANDARD DEVIATION, I.E. ALL
!	PARTICLES WERE ASSUMED TO GROW AT THE SAME RATE. THIS
!	ASSUMPTION WAS RELAXED AT VERSION 2. HOWEVER, THIS
!	INTRODUCED A SUBTLE BUG: WHEN FITZGERALD'S BETA PARAMETER
!	IS NOT UNITY, ONE HAS TO TAKE INTO ACCOUNT THE FACT THAT
!	HIS FORMULA FOR ALPHA WAS STATED IN A FORM DERIVED FOR
!	PARTICLE RADII IN MICRONS. SINCE THIS PROGRAM WORKS IN
!	METRES, ALPHA HAS TO BE INCREASED BY A FACTOR WHICH TURNS
!	OUT TO BE 10.0**( 6.0*(BETA-1) ). VERSION 3 DIFFERED FROM
!	VERSION 2 BY INCLUDING THIS FACTOR.
!       VERSION 3, HOWEVER, PRODUCED A DISCONTINUITY IN THE SIZE
!       OF THE PARTICLES AT RELATIVE HUMIDITY 0.81. VERSION 4
!       CORRECTS THIS BY INTERPOLATING BETA BETWEEN HUMIDITIES
!       0.3 AND 0.81, IN THE SAME WAY AS ALPHA HAS ALWAYS BEEN
!       INTERPOLATED. VERSION 5 IS ALMOST THE SAME AS VERSION 4,
!       BUT WITH THE UNIT CONVERSION CORRECTION ONLY APPLIED AT
!       HUMIDITIES GREATER THAN OR EQUAL TO 0.81, AS IN VERSION 3.
!
!	ONLY ONE CHEMICAL IS CATERED FOR AT PRESENT,
!	AMMONIUM SULPHATE. THUS ICHEM_TYPE IS REDUNDANT.
!
!	ABOVE THE DELIQUESCENCE POINT, TAKEN AS 0.81,
!	THE SCHEME IS THE ONE DUE TO J.W.FITZGERALD (1975).
!	IT IS NOT VERY GOOD BUT WILL DO FOR THE TIME BEING.
!	THE HYSTERESIS EFFECT IS MODELLED BY A REGION WHERE GROWTH
!	IS LINEAR. THE JUSTIFICATION FOR THIS IS THAT NEAR THE
!	DELIQUESCENCE POINT, THE PROBABILITY THAT A PARTICLE HAS
!	PREVIOUSLY BEEN EXPOSED TO A HUMIDITY IN EXCESS OF THAT
!	NEEDED FOR DELIQUESCENCE IS HIGH, SO THAT ON AVERAGE THE
!	UPPER BRANCH OF THE GROWTH CURVE SHOULD RECEIVE A HIGHER
!	WEIGHTING THAN THE LOWER BRANCH. AT LOW HUMIDITIES THE
!	CONVERSE APPLIES.
!	NOTE THAT A CEILING IS PLACED AT A HUMIDITY OF 0.995,
!	BECAUSE OF THE HIGHLY NONLINEAR GROWTH AS SATURATION
!	IS APPROACHED. IN ANY CASE, AT SUCH HIGH HUMIDITIES THE
!	IMPACT OF AEROSOL WILL BE SWAMPED BY THE EFFECT OF CLOUD.
!
!	NOTE THAT FOR THE PURPOSE OF DILUTING THE REFRACTIVE INDEX
!	OF A PARTICLE, IT IS STILL ASSUMED THAT ALL PARTICLES GROW
!	BY THE SAME FACTOR. (THIS IS RETURNED AS GROWTH_FACTOR.)
!	THIS IS NOT A SERIOUS INCONSISTENCY, BECAUSE BY THE TIME
!	BETA IS SIGNIFICANTLY DIFFERENT FROM UNITY, THE PARTICLES
!	ARE VERY DILUTE, SO THAT SMALL ERRORS IN THE AMOUNT OF WATER
!	DO NOT MAKE A LARGE DIFFERENCE TO THE REFRACTIVE INDEX.
!
!	WARNING.
!	--------
!
!	A SPECIAL PATCH OF CODE HAS BEEN ADDED TO GENERATE VALUES AT
!	A HUMIDITY OF 1.0. IT SHOULD BE NOTED THAT THIS IS OUTSIDE 
!	THE ALLOWABLE RANGE OF FITZGERALD'S SCHEME. THE VALUES 
!	RETURNED HAVE NO RELATION TO WHAT ACTUALLY HAPPENS AT 
!	SATURATION IN REALITY. THE POINT OF THIS PATCH IS THAT VALUES 
!	AT A HUMIDITY OF UNITY ARE NECESSARY FOR THE LOOK-UP TABLE.
!
!	THE VALUES OF ALPHA AND BETA AT HUMIDITY 1.0 ARE CHOSEN BY
!	EXTRAPOLATING LINEARLY FROM THE VALUES AT 0.995 AND A LOWER
!	HUMIDITY, LASTHUM, THAT SHOULD BE SET TO THE LAST HUMIDITY
!	VALUE BEFORE 1.0 IN THE LOOKUP TABLE. THE EFFECT OF THIS
!	IS THAT INTERPOLATION BETWEEN HUMIDITY VALUES LASTHUM AND 1.0
!	WILL GIVE THE VALUES OF ALPHA AND BETA THAT WOULD BE OBTAINED
!	BY LINEAR INTERPOLATION BETWEEN LASTHUM AND 0.995.
!
!	THIS IS STILL NOT EXACTLY IDEAL, BUT IT IS THE BEST I CAN DO
!	FOR THE TIME BEING.
!
!
! Current Owner of Code: D. L. Roberts/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.
!
!- ---------------------------------------------------------------------
      SUBROUTINE GROW_PARTICLES(HUMIDITY,LN_R0_LN,LN_SIGMA_LN
     &          ,ICHEM_TYPE,GROWTH_FACTOR,LASTHUM,IERR)
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
!
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/error.pcf'
      INCLUDE 'incl/aerosol_component.pcf'
!
      REAL*8
     &  HUMIDITY,
!              RELATIVE HUMIDITY
     &  LN_R0_LN,
!              LOG OF MEAN OF LOG-NORMAL DISTN.
     &  LN_SIGMA_LN,
!              LOG OF ST.DEV. OF LOG-NORMAL DISTN.
     &  GROWTH_FACTOR,
!              GROWTH FACTOR FOR DILUTION.
     &  ALPHA,
!              FITZGERALD'S ALPHA PARAMETER
     &  BETA,
!              FITZGERALD'S BETA PARAMETER
     &  PHI,
!              FITZGERALD'S PHI PARAMETER
     &  ALPHA81,
!              VALUE OF ALPHA AT HUMIDITY 0.81.
     &  BETA81,                
!              VALUE OF BETA AT HUMIDITY 0.81.
     &  POWER
!              EXPONENT FOR UNITS CONVERSION FACTOR.
!
!     DECLARATIONS FOR SPECIAL PATCH.
!
      REAL*8
     & LASTHUM,
!              LAST LOOKUP HUMIDITY VALUE BELOW 1.0
     & ALPHA_LAST,
!              VALUE OF ALPHA AT LASTHUM
     & BETA_LAST
!              VALUE OF BETA AT LASTHUM
!
      INTEGER IERR             ! ERROR FLAG
      INTEGER ICHEM_TYPE       ! INDICATOR FOR CHEMICAL TYPE.
!     ALLOWED VALUES ARE:  IP_AMMONIUM_SULPHATE.
!*******************
!these are not allowed in Ellie's version
!                          IP_ACCUM_SULPHATE.
!                          IP_AITKEN_SULPHATE.
!**************************
!
!     STATEMENT FUNCTION SPECIFICATIONS.
!
      REAL*8
     & H,
!              DUMMY VARIABLE FOR HUMIDITY
     & P,
     & AFUNC1,
     & AFUNC2,
     & BFUNC,
     & PFUNC
!
      AFUNC1(H) = 1.2D+00*EXP( (0.066D+00*H)/(1.058D+00 - H) )
      BFUNC(H) = EXP( (7.7D-04*H)/(1.009D+00 - H) )
      PFUNC(H) = 1.058D+00 - ( (0.0155D+00*(H - 0.97D+00))
     &                  / (1.02D+00 - H**1.4D+00) )
      AFUNC2(H,P) = 1.2D+00*EXP( (0.066D+00*H)/(P - H) )
!
!     SELECT TYPE OF NUCLEUS.
!
!replace original code to deal with one aerosol type only EJH 1999
!      IF ( ( ICHEM_TYPE .EQ. IP_AMMONIUM_SULPHATE ) .OR.
!     &     ( ICHEM_TYPE .EQ. IP_ACCUM_SULPHATE ) .OR.
!     &     ( ICHEM_TYPE .EQ. IP_AITKEN_SULPHATE ) ) THEN

      IF ( ( ICHEM_TYPE .EQ. IP_AMMONIUM_SULPHATE ) ) THEN
!
!     AMMONIUM SULPHATE SECTION.
!
        ALPHA = 1.0D+00
        BETA = 1.0D+00
        GROWTH_FACTOR = 1.0D+00
!
!     SPECIFY ALPHA.
!     NO GROWTH TAKES PLACE BELOW A HUMIDITY OF 0.30.
!
        IF ( HUMIDITY .GE. 0.3D+00 ) THEN
!
          IF ( HUMIDITY .LT. 0.81D+00 ) THEN
!
!     WE HAVE TO BE CAREFUL HERE TO USE THE ACTUAL VALUE OF ALPHA
!     AT 0.81, TAKING INTO ACCOUNT THE UNIT CONVERSION FACTOR.
!     OTHERWISE THE FUNCTION WOULD BECOME DISCONTINUOUS AT 0.81.
!
            POWER = 6.0D+00*( BFUNC(0.81D+00) - 1.0D+00 )
            ALPHA81 = AFUNC1(0.81D+00)*(10.0D+00**POWER)
            ALPHA = 1.0D+00 + (HUMIDITY-0.3D+00)
     &         *(ALPHA81-1.0D+00)/0.51D+00
          ELSE IF ( HUMIDITY .LE. 0.97D+00 ) THEN
            ALPHA = AFUNC1(HUMIDITY)
          ELSE IF ( HUMIDITY .LE. 0.995D+00 ) THEN
            PHI = PFUNC(HUMIDITY)
            ALPHA = AFUNC2(HUMIDITY,PHI)
          ELSE IF ( HUMIDITY .LT. 0.99999D+00 ) THEN
!
!     THE MAXIMUM VALUE OF 0.995 FOR THE HUMIDITY LEADS TO
!     THIS VALUE OF ALPHA.
!
            PHI = PFUNC(0.995D+00)
            ALPHA = AFUNC2(0.995D+00,PHI)
          ELSE
!
!     SPECIAL PATCH FOR UNIT HUMIDITY. (SEE HEALTH WARNING.)
!     NOTE THAT WE AGAIN HAVE TO PUT IN THE UNIT CONVERSION
!     FACTOR, IN ORDER TO GET THE EXTRAPOLATION RIGHT.
!
            IF ( LASTHUM .LE. 0.3D+00 ) THEN
              ALPHA_LAST = 1.0D+00
            ELSE IF ( LASTHUM .LT. 0.81D+00 ) THEN
              POWER = 6.0D+00*( BFUNC(0.81D+00) - 1.0D+00 )
              ALPHA81 = AFUNC1(0.81D+00)*(10.0D+00**POWER)
              ALPHA_LAST = 1.0D+00 + (LASTHUM-0.3D+00)
     &           *(ALPHA81-1.0D+00)/0.51D+00
            ELSE IF ( LASTHUM .LE. 0.97D+00 ) THEN
              POWER = 6.0D+00*( BFUNC(LASTHUM) - 1.0D+00 )
              ALPHA_LAST = AFUNC1(LASTHUM)*(10.0D+00**POWER)
            ELSE IF ( LASTHUM .LT. 0.995D+00 ) THEN
              PHI = PFUNC(LASTHUM)
              POWER = 6.0D+00*( BFUNC(LASTHUM) - 1.0D+00 )
              ALPHA_LAST = AFUNC2(LASTHUM,PHI)*(10.0D+00**POWER)
            ELSE
              WRITE(IU_ERR,998)
              IERR = I_ERR_FATAL
              RETURN
            ENDIF
!
!     SET ALPHA TO ITS VALUE AT 0.995 BEFORE EXTRAPOLATING TO 1.
!     NOTE THAT WE AGAIN HAVE TO PUT IN THE UNIT CONVERSION
!     FACTOR, IN ORDER TO GET THE EXTRAPOLATION RIGHT.
!
            PHI = PFUNC(0.995D+00)
            POWER = 6.0D+00*( BFUNC(0.995D+00) - 1.0D+00 )
            ALPHA = AFUNC2(0.995D+00,PHI)*(10.0D+00**POWER)
            ALPHA = ( (1.0D+00-LASTHUM)*ALPHA - 0.005D+00*ALPHA_LAST )
     &              / (0.995D+00 - LASTHUM)
          ENDIF
        ENDIF
!
!     SPECIFY BETA.
!     BELOW THE DELIQUESCENCE POINT, WE INTERPOLATE BETA
!     IN THE SAME WAY AS ALPHA. (THIS IS A CHANGE FROM
!     VERSION 3, IN WHICH IT WAS ASSUMED THAT ALL
!     PARTICLES GREW BY THE SAME FACTOR, SO BETA=1.0.)
!
        IF ( HUMIDITY .GE. 0.3D+00 .AND. HUMIDITY .LT. 0.81D+00 ) THEN
          BETA81 = BFUNC(0.81D+00)
          BETA = 1.0D+00 + (HUMIDITY-0.3D+00)*(BETA81-1.0D+00)/0.51D+00
        ENDIF
        IF ( HUMIDITY .GE. 0.81D+00 ) THEN
          BETA = BFUNC(HUMIDITY)
        ENDIF
        IF ( HUMIDITY .GT. 0.995D+00 ) THEN
          BETA = BFUNC(0.995D+00)
        ENDIF
!
!     SPECIAL PATCH FOR UNIT HUMIDITY. (SEE HEALTH WARNING.)
!
        IF ( HUMIDITY .GE. 0.99999D+00 ) THEN
          IF ( LASTHUM .LE. 0.3D+00 ) THEN
            BETA_LAST = 1.0D+00
          ELSE IF ( LASTHUM .LT. 0.81D+00 ) THEN
            BETA81 = BFUNC(0.81D+00)
            BETA_LAST = 1.0D+00 
     &        + (LASTHUM-0.3D+00)*(BETA81-1.0D+00)/0.51D+00
          ELSE IF ( LASTHUM .LT. 0.995D+00 ) THEN
            BETA_LAST = BFUNC(LASTHUM)
          ELSE
            WRITE(IU_ERR,998)
            IERR = I_ERR_FATAL
            RETURN
          ENDIF
!
!     NOW EXTRAPOLATE TO UNIT HUMIDITY.
!
          BETA = ( (1.0D+00-LASTHUM)*BETA - 0.005D+00*BETA_LAST )
     &              / (0.995D+00 - LASTHUM)
        ENDIF
!
!     INCREASE LOGS OF DISTRIBUTION MEAN AND STANDARD DEVIATION.
!
        LN_R0_LN = BETA*LN_R0_LN + LOG(ALPHA)
        LN_SIGMA_LN = BETA*LN_SIGMA_LN
        GROWTH_FACTOR = ALPHA
!
!     NOW CORRECT FOR THE UNIT CONVERSION PROBLEM WHEN HUMIDITY
!     IS AT OR ABOVE 0.81. THIS ARISES BECAUSE FITZGERALD'S
!     FORMULA IS DERIVED FOR PARTICLE RADII IN MICRONS.
!     WE DO NOT MAKE THIS CORRECTION BELOW HUMIDITY=0.81,
!     BECAUSE THE TRUE VALUES OF ALPHA (I.E. ALREADY CORRECTED)
!     WERE USED IN THE INTERPOLATION. THIS IS THE ONLY
!     SIGNIFICANT CHANGE FROM VERSION 4 TO VERSION 5.
!
        IF ( HUMIDITY.GE.0.81D+00 .AND. HUMIDITY.LT.9.9999D-01 ) THEN
          POWER = 6.0D+00*(BETA - 1.0D+00)
          LN_R0_LN = LN_R0_LN + POWER*LOG(10.0D+00)
          GROWTH_FACTOR = GROWTH_FACTOR*(10.0D+00**POWER)
        ENDIF
      ELSE
!
!     TRAP INCORRECT TYPES HERE.
!
        WRITE(IU_ERR,999) ICHEM_TYPE
        IERR = I_ERR_FATAL
!
      ENDIF
      print*,exp(LN_R0_LN),exp(LN_SIGMA_LN),GROWTH_FACTOR
!
      RETURN
!
998   FORMAT(1X,'LASTHUM should be less than 0.995')
999   FORMAT(1X,'Unavailable aerosol composition type',I5)
      END








