c     description of inputs
c     sunlat:  sun latitude, i.e., declination (degrees)
c     sunlon:  sun longitude (degrees)
c     sun_to_earth:  sun to earth distance (km)
c     hour:  GMT hour of day
c     earth_lat:  earth latitude (degrees)
c     earth_lon:  earth east longitude (degrees) 

c     description of output:
c     diurnal_sst:  diurnal warming of the sub-skin layer (K)

      subroutine get_diurnal_sst_v2_export(sunlat,sunlon,sun_to_earth,hour,earth_lat,earth_lon,wind, diurnal_sst)
      implicit none

 	real(4), parameter :: solar_constant=1361    !(joules/sec)/m**2)
 	real(4), parameter :: sun_to_earth0=149.597d6 !average distance between sun and earth
 	
	real(4), parameter :: cp=4185.5              !heat capacity joules/(kg C)
	real(4), parameter :: rho=998.               !water density kg/m**2
      real(4), parameter :: alpha_cloud=0.25       !fraction of radiation backscattered by cloud
      real(4), parameter :: alpha_gas=0.20         !Atmospheric Absorption Coefficient
	real(4), parameter :: d0=0.50357             !subskin depth meters, using exactly 0.5 changes diurnal_sst by 0.01 at wind=1 m/s, less at higher winds
	real(4), parameter :: windexp=0.46           !wind exponent term s/m
 
      real(4), parameter :: dt=0.1  !time step in hours for integration
      
	real(4) sunlat,sunlon,sun_to_earth,hour,earth_lat,earth_lon,wind
      real(4) diurnal_sst
      real(4) solar_flux_surface,coszenith
      
      integer(4) i
      real(4) sunlonx,flux,term1,term2,hourx,time_const,hour_local,dang,sunlon_sunrise
      real(4) dhour,dhr,d
      real(8) tsum
      

      term1=sind(earth_lat)*sind(sunlat)
      term2=cosd(earth_lat)*cosd(sunlat)
      
c     ================================================================================================
c     ======================================= find hour after sunrise ================================
c     ====================== this is only used for computing the early morning lag ===================
c     ================================================================================================
      if(abs(term1/term2).gt.1) then
      dhr=-999  !there is no sunrise, either all day or all  night

      else
      dang=acosd(-term1/term2)  !0 to 180
      sunlon_sunrise=earth_lon + dang !sun longitude decreases with time, thus to go back in time and find longitude at sunrise, you need to add dang
      dhr=(sunlon_sunrise-sunlon)/15.  !hour after sunrise
      if(dhr.ge.24) dhr=dhr-24.
      if(dhr.ge.24) dhr=dhr-24.
      if(dhr.lt. 0) dhr=dhr+24.
	endif
	
c     =============================================================================================
c     ================================= find cooling time constant  ===============================
c     =============================================================================================
      hour_local=hour + earth_lon/15.
      time_const=2.75
      if(hour_local.ge.7 .and. hour_local.le.19) then
      time_const=2.75
      else
      time_const=2.75 + 1.1*cosd(15.*( hour_local-1.))**2
      endif
	
	
c     =============================================================================================
c     ======================= find time indedendent part of surface flux ==========================
c     ================= factor of 3600 convert flux from power/sec to power/hour ==================
c     == atmospheric absorption varies with solar zenith angle and is applied inside the integral =
c     =============================================================================================
      solar_flux_surface=3600*(1-alpha_cloud)*solar_constant*(sun_to_earth0/sun_to_earth)**2   
     
c     =============================================================================================
c     ================================= do backward integration ===================================
c     =============================================================================================
      tsum=0
      do i=0,120
      hourx=dt*i  !hours before current time

      sunlonx=sunlon + 15.*hourx
      coszenith=term1 + term2*cosd(sunlonx-earth_lon) 
      if(coszenith.lt.0.01) cycle  !too small to worry about
      
      flux=solar_flux_surface*coszenith*exp(-alpha_gas/coszenith)
      
c     find early morning lag
      if(dhr.gt.-998) then
      dhour=dhr-hourx   !hours after sunrise
 	if(dhour.lt.0) dhour=dhour+24.
      if(dhour.le.4) flux=(1-exp(-0.25*dhour**2))*flux  !speed up program, dhour>4 represents essentially no change
      endif
     
      tsum=tsum + exp(-hourx/time_const)*flux
      enddo  !i
      
      d=d0*exp(windexp*wind)  !depth of warming layer
    
	diurnal_sst=tsum*dt/(d*rho*cp)
	
	return
	end