#include "ctrparam.h"

!	==========================================================
!
!	SULFUT.F:  Subroutine for setup prescribed sulfate
!			aerosol forcing in the format of
!			surface albedo.
!
!	----------------------------------------------------------
! 
!	Revision History:
!	
!	When	Who		What
!	-----	----------	-------	
!	080100	Chien Wang	repack based on CliChem3 & M24x11,
!				  and add cpp.
!	081100	Chien/Andrei	open cpp control for all case.
!
!	==========================================================

      subroutine sulfr_2050(BSO4LAND,BSO4OCEAN,TNOW)

#include "BD2G04.COM"
      parameter (NYS=22)

      dimension BSO4LAND(JM0),BSO4OCEAN(JM0),
     & BSO4LD86(JM0),BSO4ON86(JM0),
     & BSO4LD50(JM0),BSO4ON50(JM0),
     & IYSULF(20),CF86(20),CF50(20)
      dimension CF86N(NYS),CF50N(NYS),YSULF(NYS),
     & BSO4LD(JM0,NYS),BSO4ON(JM0,NYS)
      character * 120 sulf1986,sulf2050,sulfamp
      common /sulfdata/sulf1986,sulf2050,sulfamp
       logical first
       data first /.true./

       if(first) then
       JM=JM0
       open (575,file=sulf1986,status='old')
       open (576,file=sulf2050,status='old')
       open (577,file=sulfamp,status='old')
       read (575,'(E13.5)'),(BSO4LD86(J),j=1,JM)
       read (575,'(E13.5)'),(BSO4ON86(J),j=1,JM)
       read (576,'(E13.5)'),(BSO4LD50(J),j=1,JM)
       read (576,'(E13.5)'),(BSO4ON50(J),j=1,JM)
       read (577,'(I4)'),(IYSULF(J),j=1,20)
       read (577,'(E13.5)'),(CF86(J),j=1,20)
       read (577,'(E13.5)'),(CF50(J),j=1,20)
       close (575)
       close (576)
       close (577)
c      print *,' BSO4LD86'
c      print ('(6F7.4)'),(BSO4LD86(J)*1.e6,j=1,JM)
c      print *,' BSO4ON86'
c      print ('(6F7.4)'),(BSO4ON86(J)*1.e6,j=1,JM)
c      print *,' BSO4LD50'
c      print ('(6F7.4)'),(BSO4LD50(J)*1.e6,j=1,JM)
c      print *,' BSO4ON50'
c      print ('(6F7.4)'),(BSO4ON50(J)*1.e6,j=1,JM)
       print *,'From sulfut 2050'
       print *,' IYSULF'
       print ('(5i5)'),(IYSULF(J),j=1,20)
       print *,' CF86'
       print ('(5E13.5)'),(CF86(J),j=1,20)
       print *,' CF50'
       print ('(5E13.5)'),(CF50(J),j=1,20)
       do n=1,20
        YSULF(n)=IYSULF(n)+0.5
        do j=1,jm
         BSO4LD(j,n)=CF86(n)*BSO4LD86(J)+CF50(n)*BSO4LD50(J)
         BSO4ON(j,n)=CF86(n)*BSO4ON86(J)+CF50(n)*BSO4ON50(J)
        enddo
c      print *,' n=',n,' YSULF(n)=',YSULF(n)
c      print *,' BSO4LD'
c      print ('(6F7.4)'),(BSO4LD(J,n)*1.e6,j=1,JM)
       enddo
       if(NYS.gt.20)then
          YSULF(21)=2150.5
          YSULF(22)=2201.0
          do n=21,NYS
           do j=1,jm
             BSO4LD(j,n)=CF50(20)*BSO4LD50(J)
             BSO4ON(j,n)=CF50(20)*BSO4ON50(J)
           enddo
           print *,' n=',n,' YSULF(n)=',YSULF(n)
           print *,' BSO4LD'
           print ('(6F7.4)'),(BSO4LD(J,n)*1.e6,j=1,JM)
         enddo
       endif
       print *,' YSULF'
       print ('(5E13.5)'),(YSULF(J),j=1,NYS)
       first=.false.
       end if
       do n=1,NYS-1
        if(TNOW.gt.YSULF(n).and.TNOW.le.YSULF(n+1))go to 100
       enddo
       print *,' Wrong TNOW TNOW=',TNOW
       stop
  100  continue
        x=(YSULF(n+1)-TNOW)/(YSULF(n+1)-YSULF(n))
        do j=1,jm
         BSO4LAND(J)=x*BSO4LD(j,n)+(1.-x)*BSO4LD(j,n+1)
         BSO4OCEAN(J)=x*BSO4ON(j,n)+(1.-x)*BSO4ON(j,n+1)
        enddo
c       print *,' TNOW=',TNOW
c       print *,' n=',n,' YSULF(n)=',YSULF(n),
c    &   ' YSULF(n+1)=',YSULF(n+1)
c      print *,' BSO4LAND'
c      print ('(6F7.4)'),(BSO4LAND(J)*1.e6,j=1,JM)
c      print *,' BSO4OCEAN'
c      print ('(6F7.4)'),(BSO4OCEAN(J)*1.e6,j=1,JM)

       return
       end