C $Header: /u/gcmpack/MITgcm/pkg/dic/insol.F,v 1.14 2010/07/30 12:06:44 mlosch Exp $
C $Name:  $

#include "DIC_OPTIONS.h"

CBOP
C !ROUTINE: INSOL

C !INTERFACE: ==========================================================
      SUBROUTINE INSOL(Time,sfac,bi,bj,myThid)

C !DESCRIPTION:
C find light as function of date and latitude
C based on paltridge and parson


C !USES: ===============================================================
      IMPLICIT NONE
C     === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "FFIELDS.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "DIC_VARS.h"

C !INPUT PARAMETERS: ===================================================
C Time                 :: current time
       _RL Time
       INTEGER bi,bj
       INTEGER myThid

C !OUPUT PARAMETERS: ===================================================
       _RL sfac(1-OLy:sNy+OLy)

#ifdef DIC_BIOTIC

C !LOCAL VARIABLES: ====================================================
       _RL  solar, albedo
       _RL  dayfrac, yday, delta
       _RL  lat, sun1, dayhrs
       _RL  cosz, frac, fluxi
       integer j
CEOP

      solar = 1360. _d 0   !solar constant
      albedo= 0.6 _d 0     !planetary albedo

C     Case where a 2-d output array is needed: for now, stop here.
      IF ( usingCurvilinearGrid .OR. rotateGrid ) THEN
       STOP 'ABNORMAL END: S/R INSOL: 2-D output not implemented'
      ENDIF

C find day (****NOTE for year starting in winter*****)
        dayfrac=mod(Time,360. _d 0*86400. _d 0)
     &                    /(360. _d 0*86400. _d 0)  !fraction of year
        yday = 2. _d 0*PI*dayfrac                    !convert to radians
        delta = (0.006918 _d 0
     &         -(0.399912 _d 0*cos(yday))            !cosine zenith angle
     &         +(0.070257 _d 0*sin(yday))            !(paltridge+platt)
     &         -(0.006758 _d 0*cos(2. _d 0*yday))
     &         +(0.000907 _d 0*sin(2. _d 0*yday))
     &         -(0.002697 _d 0*cos(3. _d 0*yday))
     &         +(0.001480 _d 0*sin(3. _d 0*yday)) )
       DO j=1-OLy,sNy+OLy
C latitude in radians
          lat=YC(1,j,1,bj)*deg2rad
C     latitute in radians, backed out from coriolis parameter
C     (makes latitude independent of grid)
          IF ( usingCartesianGrid .OR. usingCylindricalGrid )
     &         lat = asin( fCori(1,j,1,bj)/(2. _d 0*omega) )
          sun1 = -sin(delta)/cos(delta) * sin(lat)/cos(lat)
          IF (sun1.LE.-0.999 _d 0) sun1=-0.999 _d 0
          IF (sun1.GE. 0.999 _d 0) sun1= 0.999 _d 0
          dayhrs = abs(acos(sun1))
          cosz = ( sin(delta)*sin(lat)+              !average zenith angle
     &            (cos(delta)*cos(lat)*sin(dayhrs)/dayhrs) )
          IF (cosz.LE.5. _d -3) cosz= 5. _d -3
          frac = dayhrs/PI                           !fraction of daylight in day
C daily average photosynthetically active solar radiation just below surface
          fluxi = solar*(1. _d 0-albedo)*cosz*frac*parfrac

C convert to sfac
          sfac(j) = MAX(1. _d -5,fluxi)
       ENDDO !j

#endif /* DIC_BIOTIC */

      RETURN
      END