clmzen.f90

CLM集合卡曼滤波数据同化算法

      SUBROUTINE clmzen (kpt, calday, dlon, dlat, cosz)

!=======================================================================
!      Source file: clmzen.f90
! Original version: Yongjiu Dai, September 15, 1999
!
! ------------------------ notes ----------------------------------
! cosine solar zenith angle from:
!    o day (1.x to 365.x), where x=0 (e.g. 213.0) denotes 00:00 at greenwich
!    o latitude,  where SH = - and NH = + 
!    o longitude, where WH = - and EH = +

! the solar declination must match that used in the atmospheric model.
! for ccm2, this code matches the ccm2 cosz to within +/- 0.0001.

! this discrepancy between clm cosz and atm cosz causes a problem.
! clm cosz may be <= zero (sun below horizon), in which case albedos
! equal zero, but atm cosz may be > zero (sun above horizon), in which
! case atm model needs albedos. no problem if atm model has sun below 
! horizon but clm has sun above horizon so long as atm solar fluxes 
! equal zero. a possible solution then is to reset points with sun 
! slightly below horizon to slightly above horizon. 

! in practice this error is not very large. e.g., if albedo error is 
! 0.0001 (atm cosz = 0.0001, clm cosz = 0) absorbed solar radiation 
! error is incident flux * 0.0001. since incident flux is itself 
! multiplied by atm cosz, incident flux is small. hence, error is small.
! in fact the error is smaller than the difference between atm net solar 
! radiation at surface and clm net solar radiation at surface, which arises
! due to the different surface radiation parameterizations.

! however, reset points as discussed above just in case atm model may
! blow up if albedos equal zero when atm cosz > 0.

! Compute solar distance factor and cosine solar zenith angle usi
! day value where a round day (such as 213.0) refers to 0z at
! Greenwich longitude.
!
! Use formulas from Paltridge, G.W. and C.M.R. Platt 1976: Radiative
! Processes in Meterology and Climatology, Elsevier Scientific
! Publishing Company, New York  p. 57, p. 62,63.

! -----------------------------------------------------------------

      IMPLICIT NONE

! ------------------------ input/output variables -----------------
! input
      integer, INTENT(in) :: &
            kpt             !number of land points, including subgrid

      real, INTENT(in) :: &
            calday,        &!calendar day + fraction (1.xx -> 365.xx)
            dlat(kpt),     &!latitude  (radians): + = NH
            dlon(kpt)       !longitude (radians): + = EH
! output
      real, INTENT(out) :: &
            cosz(kpt)       !cosine zenith angle 

! ------------------------ local variables ------------------------
      real dayspy           !days per year
      real pie              !pie
      real theta            !earth orbit seasonal angle in radians
      real delta            !solar declination angle  in radians
      real sind             !sine   of declination
      real cosd             !cosine of declination
      real phi              !greenwich calendar day + longitude offset
      real loctim           !local time (hour)
      real hrang            !solar hour angle, 24 hour periodicity (radians)
      real mcsec            !current seconds in day (0, ..., 86400)

      integer k             !index
! -----------------------------------------------------------------

      dayspy = 365.
      pie = 4.*atan(1.)

! solar declination: match ccm2  

      theta = (2.*pie*calday)/dayspy
      delta = .006918 - .399912*cos(   theta) + .070257*sin(   theta) &
                      - .006758*cos(2.*theta) + .000907*sin(2.*theta) &
                      - .002697*cos(3.*theta) + .001480*sin(3.*theta)
      sind = sin(delta)
      cosd = cos(delta)

      do k = 1, kpt

! local time

         mcsec = (calday - int(calday)) * 86400.
         phi = calday + dlon(k)/(2.*pie)
         loctim = (mcsec + (phi-calday)*86400.) / 3600.
         if(loctim>=24.) loctim = loctim - 24.

! hour angle

         hrang = 15. * (loctim-12.) * pie/180.     ! 360/24 = 15

! cosine solar zenith angle. reset points with sun slightly below horizon 
! to slightly above horizon, as discussed in notes.

         cosz(k) = sin(dlat(k))*sind + cos(dlat(k))*cosd*cos(hrang)
         if (cosz(k) .ge. -0.001 .and. cosz(k).le. 0.) cosz(k) = 0.001

      end do

      END SUBROUTINE clmzen