SUBROUTINE SUN_RADFIX(SIGMA,COSAS1,COSAS2,PSOL,HEAT,TOIT) 2

*     AUTHOR:    FRANCOIS LEMAY (CMDN), SEPTEMBRE 2003
*
*     REVISIONS :
*
*     001        A-M.LEDUC (Sept. 2003)   - IBM conversion
*
*     OBJECT: THE SOLAR HEATING IS FIXED FOR THE STRATOSPHERIC LEVELS
*             (ABOVE 85MB). A NEW ROUTINE IS CREATED FOR THIS PURPOSE
*             RATHER THAN HAVING THE SAME RADFIX EXPRESSION IN
*             BOTH NEWRAD (non-radiative timesteps) AND SUN (radiative
*             timesteps) ROUTINES.

#include "impnone.cdk"

      REAL    SIGMA,SIGMA2,COSAS1,COSAS2,PSOL,P,HEAT,VV,VV2,A,ZZ,Y
      LOGICAL TOIT
      REAL    REC_86400

*     SIGMA:     sigma coordinate
*     COSAS1,2:  cosine of the solar zenith angle
*     PSOL:      surface pressure [N/m2]
*     HEAT:      heating rate [Kelvin/second]
*     TOIT:      .true. if the sigma value corresponds to the roof of the model

*     Note: the reason there is two cosines of the sun angle comes from
*           the newrad3.ftn code where the average of the cosine between
*           two radiative timesteps (cosas1) and the cosine itself (cosas2)
*           are involved in the radfix code. To get bit validation with the old code,
*           we introduced two variables but this is not very relevant.

*     LES TAUX SONT FIXES ET NON CALCULES AU DESSUS DE 85 MB

      VV=AMAX1(SIGMA,0.002)
*     vertical coordinate independant of the surface pressure: 
      VV2=VV*PSOL*1.E-5
*     VV2=VV

*     EXPOSANT VARIANT AVEC HAUTEUR AU-DESSUS DE 85MB
      
      A=0.4

*     vertical coordinate independant of the surface pressure:
      SIGMA2=SIGMA*PSOL*1.E-5
*     SIGMA2=SIGMA

*     A=CVMGT(0.6,A,SIGMA2.LT.0.025)
*     A=CVMGT(.4+.274*(0.02-SIGMA2)/.02,A,SIGMA2.LT.0.020)
      IF (SIGMA2.LT.0.020) A=.4+.274*(0.02-SIGMA2)*50
*     ZZ=CVMGT(COSAS1**A,COSAS1,SIGMA2.LT.0.085.AND.SIGMA2.GT.0.0019)
      IF (SIGMA2.LT.0.085.AND.SIGMA2.GT.0.0019) THEN 
*        ZZ=COSAS1**A
         ZZ=exp(A*log(COSAS1))
      ELSE
         ZZ=COSAS1
      ENDIF

      REC_86400=1./86400.

*     LE FIT EST DU SECOND DEGRE EN 1/SIGMA POUR SOLEIL AU NADIR
*     Y= (0.327+0.0642/VV -7.082E-5/VV**2)/86400. *ZZ

*     THE NEW FIT IS A LINEAR FUNCTION WITH A FLOOR VALUE EXPRESSED
*     IN TERMS OF PRESSURE RATHER THAN SIGMA (TOPOGR. INDEPENDANT)
*     Y= MAX((-75.*VV  + 5.25),2.0)/86400. *ZZ
*     Y= MAX((-75.*VV2 + 5.25),2.0)/86400. *ZZ
      Y= MAX((-75.*VV2 + 5.25),2.0)*REC_86400 *ZZ

*     SOLAR HEATING IS FIXED AT THE TOP
*     IF (TOIT) Y= (4.5)/86400. *ZZ
      IF (TOIT) Y= (4.5)*REC_86400 *ZZ

*     Y=CVMGT(Y,0.,COSAS2.GT.0.00001)
      IF (COSAS2.LE.0.00001) Y=0.
*     HEAT=CVMGT(Y,HEAT,SIGMA2.LT.0.085)
      IF (SIGMA2.LT.0.085) HEAT=Y

      HEAT=AMAX1(HEAT,0.)

      RETURN
      END