! Land surface model with atmosphere but only for heat
! 2006.Oct  I. Tamagawa
!
!	first version 2006.10.30
!         read atmospheric data from 'atmosdata.csv'
!					( see : subroutine getatmos )
!	  atmosdata.csv is assumed as below
!           year month day hour minute WindSpeed AirTemperature RH Precipitation SolorRadiation
!
!
!         soil parameters are set in subroutine setinit
!	  surface parameters are set in subroutine setsurface
!         initial values are set setinit
!    ABOVE 3 SUBROUTINE MUST BE CHANGED FOR YOUR SIMULATION
!
!	  Major variables are defined in typedef module, most of 
!	   subroutines and main program includes it by use typedef
!

module  typedefine
implicit none
  type soilparam
   real lambda,  spec_heatcap, rho, dz, z, thickness
  end type soilparam

  type atmosparam
   real temp, windsp, rh, qv,  qsd, qld, rho, press, prec
   real time
  end type atmosparam

  type surfaceparam
   real albedo, beta, z0m, z0t, zt, zu
  end type surfaceparam

  type work
    real a, b, c, d
  end type work

  type variables
    real temp
  end type variables

  real, parameter     :: Cp = 1004.0
  real, parameter     :: Latent = 2.5e6
  real, parameter     :: StefanBoltzmann  = 5.6704e-8
  real, parameter     :: Grav  = 9.8
  real, parameter     :: Karman  = 0.4
  real, parameter     :: Rd = 287.1

  character :: atmosfile*1024

  type datetime
   integer year, month, day, hour, minute, sec
   integer julianday
  end type datetime

  integer, parameter  :: atmosunit=11

end module typedefine



!
!   Main program
!
!

program main
  use typedefine

  real :: dt=10.0		! unit : sec
  integer, parameter ::  n=100, iend=360*264  ! unit : number


  type(variables), dimension(n) ::  ground
  type(soilparam), dimension(n) :: soil
  type(atmosparam)              :: atmos
  type(surfaceparam)            :: surf
  type(work), dimension(n)      :: tmp                    
  type(datetime)      :: start, now


  real :: t, dmpt=7200.0
  integer :: i


  write(0, *) 'Input start date and time'
  write(0, *) 'year month day hour minute second'
  read(*,*) start%year, start%month, start%day, start%hour, start%minute, start%sec
  write(0, *) start%year, start%month, start%day, start%hour, start%minute, start%sec



  call date2t( start, start, t )
  call setsurface( surf, start )
  call setinit( soil, ground, n, start )

  call  output( t, dmpt, ground, soil, n, start )
  
  do i=1, iend
    write(0,*) 'Time=', t
   call onestep( t, dt, ground, soil, atmos, surf, n, tmp, start )
   t = dt*i
    call  output( t, dmpt, ground, soil, n, start )
end do
  
  stop
end



subroutine setsurface( surf, start )
  use typedefine
  type(surfaceparam)  :: surf
  type(datetime)      :: start


  surf%albedo = 0.191
  surf%beta   = BETA    !  should be set some value
  surf%z0m    = 0.001
  surf%z0t    = 0.01

  surf%zt     = 1.5 	! Height of temp
  surf%zu     = 1.5	! Height of Wind

  return
end


subroutine setinit( soil, ground, n, start)
  use typedefine

  integer n
  type(variables), dimension(n) ::  ground
  type(soilparam), dimension(n) ::  soil
  type(datetime)      :: start
  
  integer i
  real l1i, l2i, l3i

! Soil properties, should be set
  do i=1, n
    soil(i)%thickness = 0.005
  end do

! Set z from thickness

  soil(1)%z = (soil(1)%thickness)/2.0
  do i=2, n
   soil(i)%z  = soil(i-1)%z + (soil(i)%thickness)/2.0 &
   					+ (soil(i-1)%thickness)/2.0
end do
  soil(1)%dz = soil(1)%z - 0.0

  do i=2, n
   soil(i)%dz  = soil(i)%z - soil(i-1)%z
end do

! Set temp

  do i=1, n
    ground(i)%temp = 313.49
 end do
  do i=1, n
   if (soil(i)%z < 0.05) then
    soil(i)%lambda  = 1.037
    soil(i)%spec_heatcap = 928
    soil(i)%rho     = 1897
   ELSE if ((soil(i)%z > 0.05).AND.(soil(i)%z < 0.1)) then
    soil(i)%lambda = 1.593
    soil(i)%spec_heatcap = 782
    soil(i)%rho = 2384
   ELSE
    soil(i)%lambda = 1.59
    soil(i)%spec_heatcap = 1359
    soil(i)%rho = 1748 
   end if
end do

  return
end

subroutine  output( t, dmpt, ground, soil, n, start )
  use typedefine

  integer n
  real                          :: t, dmpt
  type(variables), dimension(n) ::  ground
  type(soilparam), dimension(n) :: soil
  type(datetime)      :: start, now
  real, save                    :: next=0
  
  integer  ::  i, jj, kk

  if( t >= next ) then
      do i=1,n
        write(*,*) t,  soil(i)%z, ground(i)%temp
    end do
   next = next + dmpt
  endif

  return
end
  

subroutine onestep( t, dt, ground, soil, atmos, surf, n, tmp, start )
  use typedefine

  integer  ::  n
  type(variables), dimension(n) ::  ground
  real  ::  t, dt
  type(soilparam), dimension(n)  :: soil
  type(atmosparam)               :: atmos
  type(surfaceparam)             :: surf
  type(work), dimension(n)       :: tmp
  type(datetime)      :: start, now
  
  integer ::  i
  real    :: h, le, g
  real    :: l1i, l2i, l3i, lm, lm1, fact


  
! boundary condition 

  call getatmos( t, atmos, start)
!   write(0,*) 'ATMOS', atmos%time, atmos%temp, atmos%press, atmos%windsp
!   write(0,*) '0 : Ts = ', ground(1)%temp
  call get_hle( h, le, ground(1)%temp, atmos, surf )

  g = (1.0 - surf%albedo) * atmos%qsd + atmos%qld  &
                            - StefanBoltzmann*(ground(1)%temp**4) -h -le

!   write(0,*) 'Qsd=', atmos%qsd, "Qld=", atmos%qld
!    write(0,*) 'H = ', h, 'LE=', le, 'G =', g, 'Qnet = ', &
! 	(1.0 - surf%albedo) * atmos%qsd + atmos%qld- StefanBoltzmann*(ground(1)%temp**4) &
!      , 'Qsd', atmos%qsd, 'Qld', atmos%qld, 'Qlu', StefanBoltzmann*(ground(1)%temp**4) 

! Set matrix  a * t(i-1) + b * t(i) + c * t(i+1) = d 

! Top

  i=1
  fact = dt/(soil(i)%rho*soil(i)%spec_heatcap*soil(i)%thickness)
  l1i = soil(i)%thickness/(2.0*soil(i)%lambda)
  l2i = soil(i+1)%thickness/(2.0*soil(i+1)%lambda)
  l3i = l1i+l2i
  lm = -(soil(i)%thickness+soil(i+1)%thickness)/l3i/soil(i)%dz


   tmp(i)%b = 1.0  - fact*lm
   tmp(i)%c = fact  *  lm
   tmp(i)%d = ground(i)%temp + fact * g

  do i=2,n-1
! mean heat diffusion coeff

  lm1 = lm
  l1i = soil(i)%thickness/(2.0*soil(i)%lambda)
  l2i = soil(i+1)%thickness/(2.0*soil(i+1)%lambda)
  l3i = l1i+l2i
  lm = -(soil(i)%thickness+soil(i+1)%thickness)/l3i/soil(i)%dz
  fact = dt/(soil(i)%rho*soil(i)%spec_heatcap*soil(i)%thickness)
  
   tmp(i)%a = fact  *  lm1
   tmp(i)%b = 1.0  -fact * lm1 - fact*lm
   tmp(i)%c = fact  *  lm
   tmp(i)%d = ground(i)%temp

  end do

! Bottom ( no flux condition )
 i = n
  lm1 = lm
  fact = dt/(soil(i)%rho*soil(i)%spec_heatcap*soil(i)%thickness)
  tmp(i)%a = fact  *  lm1
  tmp(i)%b = 1.0  -fact * lm1
  tmp(i)%d = ground(i)%temp


! Tomas algorism

  do i=1,n-1
   fact =  -tmp(i+1)%a/tmp(i)%b
!  tmp(i+1)%a = tmp(i+1)%a + fact * tmp(i)%b
   tmp(i+1)%b = tmp(i+1)%b + fact * tmp(i)%c
   tmp(i+1)%d = tmp(i+1)%d + fact * tmp(i)%d
  end do

  i=n
    ground(i)%temp = tmp(i)%d/tmp(i)%b
  do i=n-1, 1, -1
    ground(i)%temp = (tmp(i)%d-tmp(i)%c*ground(i+1)%temp)/tmp(i)%b
  end do

  
  return
end

subroutine getatmos( t, atmos, start )
  use typedefine
  type(atmosparam)               :: atmos
  real t
  type(datetime)      :: start, now

  type(atmosparam), save         :: aprev, anext


  integer, save  ::  open=0

! Atmospheric data are read from atmosdata.csv

  if( open==0 ) then
    open(unit=atmosunit, file='atmosdata.csv')
    write(0,*)'open atmospheric data'
    call readoneline( aprev, start )
    call readoneline( anext, start )
    open = 1
  end if


 do while( t > anext%time ) 
   call copyatmos( anext, aprev )
   call readoneline( anext, start )
 end do

 call interporate( t, atmos, aprev, anext )

end

subroutine copyatmos( anext, aprev )
  use typedefine
  type(atmosparam)               :: aprev, anext

  aprev = anext
end


subroutine interporate( t, atmos, aprev, anext )
  use typedefine
  type(atmosparam)               :: atmos, aprev, anext
  real                           :: t

  real :: m, n, l

   l = anext%time - aprev%time
   m =  (anext%time -t)/l
   n =  (t - aprev%time)/l

   atmos%time = t
   atmos%temp = m*aprev%temp + n*anext%temp
   atmos%press = m*aprev%press + n*anext%press
   atmos%windsp = m*aprev%windsp + n*anext%windsp
   atmos%qsd = m*aprev%qsd + n*anext%qsd
   atmos%qld = m*aprev%qld + n*anext%qld
   atmos%rh = m*aprev%rh + n*anext%rh
   atmos%qv = m*aprev%qv + n*anext%qv
   atmos%rho = m*aprev%rho + n*anext%rho
   atmos%prec = m*aprev%prec + n*anext%prec


end


subroutine readoneline(anext, start )
  use typedefine
  type(atmosparam)               :: anext
  type(datetime)      :: start

  type(datetime)      :: now
  real downwardlongwave



  read(atmosunit, *) now%year, now%month, now%day, &
  			now%hour, now%minute, now%sec,      &
  	anext%windsp, anext%temp, anext%rh, anext%qsd 

!   write(0, *) now%year, now%month, now%day, &
!   			now%hour, now%minute,      &
!   	anext%windsp, anext%temp, anext%rh, anext%prec, anext%qsd 
!
!   write(0,*) 'Temp=', anext%temp, 'RH=', anext%rh
!  now%sec = 0
   
  call date2t( start, now, anext%time )

  anext%press = 101300.0
  anext%temp = anext%temp + 273.15			! C to K
  anext%rh   = 0.01*anext%rh				! 0-100 --> 0-1
  anext%qld  = downwardlongwave( anext%temp, anext%rh )
  anext%qsd  = 1.0e6 * anext%qsd / 3600.0
  anext%qv   = anext%rh * qsat( anext%temp, anext%press )
  anext%rho = (1.0/(1.0-anext%qv)) * (anext%press )/(Rd*anext%temp)

!  write(0,*) 'READONELINE Qld=', anext%qld, 'Qsd=', anext%qsd, 'Temp', anext%temp

  return
end

subroutine get_hle( h, le, ts, atmos, surf )
  use typedefine
  type(atmosparam)               :: atmos
  type(surfaceparam)             :: surf
  real h, le, ts

  real cd, ch

  call bulkcoeff( cd, ch, atmos, ts, surf )

  h = atmos%rho * Cp * ch * atmos%windsp * ( ts - atmos%temp )
  le = atmos%rho * Latent * ch * surf%beta * atmos%windsp * &
	  ( qsat( ts, atmos%press ) - atmos%qv )

!   write(0,*) 'H=', h, 'le=', le
  return
end


subroutine bulkcoeff( cd, ch, atmos, ts, surf )
  use typedefine
  type(atmosparam)               :: atmos
  type(surfaceparam)             :: surf
  real                           :: cd, ch, ts

  real                           :: ri		! Richardson number
  real                           :: tm, dudz, c1, c2, fm, fh
  real, parameter                :: b=5, c=5, d=5

  tm = (atmos%temp + ts)/2
  if( atmos%windsp > 0.2 ) then
	  dudz = atmos%windsp/surf%zu
  else
	dudz = 0.2/surf%zu
  endif

  ri = Grav * ((atmos%temp - ts)/surf%zt) / (tm * dudz*dudz)

!   write(0,*) ' Ri ', ri, ' Ta', atmos%temp, ' Ts', ts

  c1 = Karman/log( surf%zu / surf%z0m )
  c2 = Karman/log( surf%zt / surf%z0t )
  
  if(ri<0) then
    fm = 1.0 - 2*b*ri/(1+3*b*c*c1*c1*sqrt(-ri*(surf%zu)/(surf%z0m)))
    fh = 1.0 - 3*b*ri/(1+3*b*c*c1*c2*sqrt(-ri*(surf%zt)/(surf%z0t)))
  else 
    fm = 1.0/(1+2*b*ri*sqrt(1.0+d*ri))
    fh = 1.0/(1+3*b*ri*sqrt(1.0+d*ri))
  end if

   cd = c1 * c1 * fm
   ch = c1 * c2 * fh

!   write(*,*) ' Cd ', cd, ' Ch', ch

  return
end


real function qsat( temp, press )
 real temp, press
 real, parameter  ::  eps = 0.622
 real e, tc, rhov, rhod, ep


 e = esat( temp )
  ep = e/press
  qsat = eps * ep / ( 1 - (1-eps)*ep )
  return
end

real function esat( temp )
 real :: temp
 real :: tc
 real, parameter  ::  t_ice = 273.15

 tc = temp - t_ice
 if( tc >= 0  ) then
   esat = 610.78 * (10.0**(7.5*tc/(237.3+tc)))
 else
   esat = 610.78 * (10.0**(9.5*tc/(265.3+tc)))
 endif
 return
end

real function downwardlongwave( temp, rh ) 
  use typedefine

  real temp, rh, ea

  integer k
  real cdall, t4, cdp, esat
  real a, b, c

 real, dimension(21) :: cdl=(/0.0, 0.0, 0.0, 0.0, 0.0, &
                              0.0, 0.0, 0.0, 0.0, 0.0, &
			      0.0, 0.0, 0.07, 0.11, 0.19, &
			      0.4, 0.85, 1.0, 1.0, 1.0,  1.0/)
 real, dimension(21) :: cdm=(/0.0, 0.0, 0.0, 0.0, 0.0, &
                              0.0, 0.0, 0.0, 0.0, 0.0, &
			      0.0, 0.05, 0.12, 0.30, 0.40,&
			      0.50, 0.70, 0.95, 1.0, 1.0, 1.0/)
 real, dimension(21) :: cdh=(/0.0, 0.0, 0.0, 0.0, 0.0,&
                              0.0, 0.0, 0.0, 0.0, 0.05, &
			      0.15, 0.30, 0.55, 0.75, 0.87,&
			      0.95, 1.00, 1.00, 1.0, 1.0,  1.0/)
 
 

! Cloud amount by Ohno & Isa 1984 same as CReSS


  k = nint(1.0+20.0*rh)
  cdall = cdl(k)+cdm(k)+cdh(k)
  
  t4 = temp*temp
  t4 = t4*t4

  ea = rh * esat( temp );


 if( cdall > 0 ) then
    cdp = (cdl(k)+0.85*cdm(k)+0.5*cdh(k))
 else
    cdp = 0
 endif

  a = -0.49 + 0.0066*sqrt(ea)
  b = 1-(0.75-0.00005*ea)* cdp

  downwardlongwave = 0.95*StefanBoltzmann*t4 * & 
  	(1.0 + a*b)
!  write(*,*) 'T=', temp, 'RH=', rh, 'Qld=', downwardlongwave, 'Dry (0%) Qld=',  &
!		0.95*StefanBoltzmann*t4 *  (1.0 + a)
!  write(*,*) 'RH=', rh, 'Qld', a, b, 0.95*StefanBoltzmann*t4, downwardlongwave
  return 
end




integer function julday(iyyy, mm, id) 
! From numerical recipe
  implicit none
  integer, intent(in) :: mm,id,iyyy
  integer, parameter  :: igreg=15+31*(10+12*1582)
  integer             :: ja,jm,jy

   jy=iyyy
   if (jy == 0) then
   	write(0,*) ' Error no 0 year '
	stop
   end if
   if (jy < 0) jy=jy+1
   if (mm > 2) then
     jm=mm+1 
    else 
     jy=jy-1 
     jm=mm+13 
   end if 
   julday=floor(365.25*jy)+floor(30.6001*jm)+id+1720995 
   if (id+31*(mm+12*iyyy) >= igreg) then
     ja=floor(0.01*jy) 
     julday=julday+2-ja+floor(0.25*ja) 
 end if 
end function julday 

subroutine caldat(julian,mm,id,iyyy) 
! From numerical recipe
  implicit none 
  integer, intent(in) :: julian 
  integer, intent(out) :: mm,id,iyyy 


  integer :: ja,jalpha,jb,jc,jd,je 
  integer, PARAMETER :: igreg=2299161 
  if (julian >= igreg) then
     jalpha=int(((julian-1867216)-0.25)/36524.25)
     ja=julian+1+jalpha-int(0.25*jalpha) 
  else if (julian < 0) then
     ja=julian+36525*(1-julian/36525)
  else
     ja=julian 
  end if 

  jb=ja+1524 
  jc=int(6680.0+((jb-2439870)-122.1)/365.25) 
  jd=365*jc+int(0.25*jc) 
  je=int((jb-jd)/30.6001) 
  id=jb-jd-int(30.6001*je) 
  mm=je-1 
  if (mm > 12) mm=mm-12 
  iyyy=jc-4715 
  if (mm > 2) iyyy=iyyy-1 
  if (iyyy <= 0) iyyy=iyyy-1 
  if (julian < 0) iyyy=iyyy-100*(1-julian/36525) 
end subroutine caldat 



subroutine setnow( t, start, now )
! in    t (sec) and start
! set    now

  use typedefine
  real               :: t
  type(datetime)     :: start, now
  integer, parameter ::  dayinsec=86400
  integer t1, jadd

  t1 = 0
  if(t>dayinsec) then
   	t1 = t/dayinsec
	t = t - t1*dayinsec
  end if

  now%julianday = start%julianday + t1

  call caldat(now%julianday,now%month, now%day, now%year) 
 
  t1 = 0
  if(t>3600) then
       t1 = t/3600
       t  = t-t1*3600
  endif 
  now%hour = t1

  t1 = 0
  if(t>60) then
       t1 = t/60
       t  = t-t1*60
  endif 
  now%minute = t1

  now%sec = t

end




subroutine date2t( start, now, t )
  use typedefine
  type(datetime) start, now
  real t

  integer jj
  integer julday
  integer, parameter ::  dayinsec=86400

  if(start%julianday == 0 ) then 
     start%julianday = julday(start%year, start%month, start%day)
  end if

     now%julianday = julday(now%year, now%month, now%day)

 t = dayinsec * (now%julianday-start%julianday)  &
    + 3600.0 * (now%hour - start%hour)		&
    + 60.0   * (now%minute - start%minute)	&
    +          ( now%sec - start%sec )

end 
