pkg/exf/exf_zenithangle.F

94390f4f16 Gael*0001 #include "EXF_OPTIONS.h"
                0002
9675e9700b Jean*0003 CBOP
                0004 C     !ROUTINE: EXF_ZENITHANGLE
                0005 C     !INTERFACE:
                0006       SUBROUTINE EXF_ZENITHANGLE(
                0007      I               myTime, myIter, myThid )
                0008
                0009 C     !DESCRIPTION: \bv
                0010 C     *==========================================================*
                0011 C     | SUBROUTINE EXF_ZENITHANGLE
                0012 C     | o compute zenith angle, derive albedo and
                0013 C     |   the incoming flux at the top of the atm.
                0014 C     *==========================================================*
                0015 C     \ev
                0016
                0017 C     !USES:
94390f4f16 Gael*0018       IMPLICIT NONE
                0019
                0020 C     == global variables ==
                0021 #include "EEPARAMS.h"
                0022 #include "SIZE.h"
                0023 #include "PARAMS.h"
                0024 #include "GRID.h"
                0025 #include "EXF_PARAM.h"
                0026 #include "EXF_FIELDS.h"
                0027 #include "EXF_CONSTANTS.h"
9675e9700b Jean*0028 #ifdef ALLOW_CAL
94390f4f16 Gael*0029 # include "cal.h"
9675e9700b Jean*0030 #endif
94390f4f16 Gael*0031
9675e9700b Jean*0032 C     !INPUT/OUTPUT PARAMETERS:
94390f4f16 Gael*0033       _RL     myTime
                0034       INTEGER myIter
                0035       INTEGER myThid
                0036
                0037 #ifdef ALLOW_DOWNWARD_RADIATION
                0038 #ifdef ALLOW_ZENITHANGLE
9675e9700b Jean*0039 C     !FUNCTIONS:
                0040 #ifdef ALLOW_CAL
                0041       INTEGER  cal_IsLeap
                0042       EXTERNAL cal_IsLeap
                0043 #endif
                0044
                0045 C     !LOCAL VARIABLES:
                0046       INTEGER i, j, bi, bj
a303c567c2 Jean*0047       INTEGER iLat1,iLat2,iTyear1,iTyear2
e19d0cc6b2 Gael*0048       _RL     wLat1,wLat2,wTyear1,wTyear2
94390f4f16 Gael*0049       _RL H0, dD0dDsq, CZENdaily, CZENdiurnal
                0050       _RL TDAY, TYEAR, ALBSEA1, ALPHA, CZEN, CZEN2
                0051       _RL DECLI, ZS, ZC, SJ, CJ, TMPA, TMPB, TMPL, hlim
9675e9700b Jean*0052       _RL SOLC, FSOL
                0053 c     _RL CSR1, CSR2, FLAT2
                0054       _RL secondsInYear
                0055 #ifdef ALLOW_CAL
                0056       _RL myDateSeconds
a303c567c2 Jean*0057       INTEGER year0,mydate(4),difftime(4)
                0058       INTEGER dayStartDate(4),yearStartDate(4)
9675e9700b Jean*0059 #endif
                0060 CEOP
94390f4f16 Gael*0061
a303c567c2 Jean*0062 C solar constant
                0063 C --------------
94390f4f16 Gael*0064       SOLC   = 1368. _d 0
a303c567c2 Jean*0065 C note: it is fourth (342. _d 0) is called SOLC in pkg/aim_v23
94390f4f16 Gael*0066
a303c567c2 Jean*0067 C determine time of year/day
                0068 C --------------------------
94390f4f16 Gael*0069
9675e9700b Jean*0070 #ifdef ALLOW_CAL
                0071       IF ( useCAL ) THEN
4d1cab9739 Jean*0072        CALL cal_GetDate( myIter, myTime, mydate, myThid )
                0073        year0         = INT(mydate(1)/10000.)
94390f4f16 Gael*0074        secondsInYear = ndaysnoleap * secondsperday
a303c567c2 Jean*0075        IF ( cal_IsLeap(year0,myThid) .EQ. 2)
94390f4f16 Gael*0076      &       secondsInYear = ndaysleap * secondsperday
a303c567c2 Jean*0077
94390f4f16 Gael*0078        yearStartDate(1) = year0 * 10000 + 101
                0079        yearStartDate(2) = 0
                0080        yearStartDate(3) = mydate(3)
                0081        yearStartDate(4) = mydate(4)
                0082        CALL cal_TimePassed(yearStartDate,mydate,difftime,myThid)
                0083        CALL cal_ToSeconds (difftime,myDateSeconds,myThid)
a303c567c2 Jean*0084
94390f4f16 Gael*0085        TYEAR=myDateSeconds/secondsInYear
a303c567c2 Jean*0086
94390f4f16 Gael*0087        dayStartDate(1) = mydate(1)
                0088        dayStartDate(2) = 0
                0089        dayStartDate(3) = mydate(3)
                0090        dayStartDate(4) = mydate(4)
                0091        CALL cal_TimePassed(dayStartDate,mydate,difftime,myThid)
                0092        CALL cal_ToSeconds (difftime,myDateSeconds,myThid)
a303c567c2 Jean*0093
94390f4f16 Gael*0094        TDAY= myDateSeconds / ( 86400 . _d 0 )
a303c567c2 Jean*0095
9675e9700b Jean*0096       ELSE
                0097 #else /* ALLOW_CAL */
                0098       IF ( .TRUE. ) THEN
                0099 #endif /* ALLOW_CAL */
                0100
                0101        secondsInYear = 86400. _d 0 * 365.25 _d 0
                0102        TYEAR = myTime / secondsInYear
                0103        TDAY  = myTime / 86400 . _d 0
                0104        TYEAR = MOD( TYEAR, oneRL )
                0105        TDAY  = MOD( TDAY , oneRL )
94390f4f16 Gael*0106
9675e9700b Jean*0107       ENDIF
                0108
                0109       IF ( useExfZenAlbedo ) THEN
d106b5e2d8 Gael*0110
94390f4f16 Gael*0111        DO bj = myByLo(myThid),myByHi(myThid)
                0112         DO bi = myBxLo(myThid),myBxHi(myThid)
                0113
9675e9700b Jean*0114          IF ( select_ZenAlbedo.EQ.0 ) THEN
f71b2ceff6 Gael*0115
9675e9700b Jean*0116           DO j = 1,sNy
                0117            DO i = 1,sNx
                0118             zen_albedo (i,j,bi,bj) = exf_albedo
                0119            ENDDO
                0120           ENDDO
f71b2ceff6 Gael*0121
9675e9700b Jean*0122          ELSEIF ( select_ZenAlbedo.EQ.1 ) then
94390f4f16 Gael*0123
a303c567c2 Jean*0124 C     This is the default option: daily mean albedo (i.e. without diurnal
                0125 C     cycle) obtained from the reference table that was computed in
                0126 C     exf_zenithangle_table.F. Using either daily or 6 hourly fields, this
                0127 C     option yields correct values of daily upward sw flux.
                0128 C     This is not the case for select_ZenAlbedo.GT.1 (see comments below).
94390f4f16 Gael*0129
                0130           iTyear1= 1 + 365.*TYEAR
                0131           wTyear1= iTyear1 - 365.*TYEAR
e19d0cc6b2 Gael*0132           iTyear2= iTyear1 + 1
                0133           wTyear2= 1.0 _d 0 - wTyear1
a303c567c2 Jean*0134
9675e9700b Jean*0135           DO j = 1,sNy
                0136            DO i = 1,sNx
                0137             IF ( zen_albedo_pointer(i,j,bi,bj).EQ. 181. _d 0 ) THEN
                0138               iLat1=181
                0139               wLat1=0.5  _d 0
                0140               iLat2=181
                0141               wLat2=0.5  _d 0
                0142             ELSE
eabfe867d8 Gael*0143               iLat1= INT(zen_albedo_pointer(i,j,bi,bj))
9675e9700b Jean*0144               wLat1= 1. _d 0 + iLat1 - zen_albedo_pointer(i,j,bi,bj)
                0145               iLat2= iLat1 + 1
                0146               wLat2= 1. _d 0 - wLat1
                0147             ENDIF
                0148             ALBSEA1 =
                0149      &         wTyear1*wLat1*zen_albedo_table(iTyear1,iLat1)
                0150      &       + wTyear1*wLat2*zen_albedo_table(iTyear1,iLat2)
                0151      &       + wTyear2*wLat1*zen_albedo_table(iTyear2,iLat1)
                0152      &       + wTyear2*wLat2*zen_albedo_table(iTyear2,iLat2)
                0153
                0154 C determine overall albedo: approximation: half direct and half diffus
                0155             zen_albedo (i,j,bi,bj) =
                0156      &          0.5 _d 0 * exf_albedo + 0.5 _d 0 * ALBSEA1
                0157            ENDDO
                0158           ENDDO
                0159
                0160 C        if select_ZenAlbedo = 0, = 1, else
                0161          ELSE
94390f4f16 Gael*0162
a303c567c2 Jean*0163 C determine solar declination
                0164 C ---------------------------
                0165 C       (formula from Hartmann textbook, after Spencer 1971)
9675e9700b Jean*0166           ALPHA= 2. _d 0*PI*TYEAR
                0167           DECLI = 0.006918 _d 0
                0168      &       - 0.399912 _d 0 * COS ( 1. _d 0 * ALPHA )
                0169      &       + 0.070257 _d 0 * SIN ( 1. _d 0 * ALPHA )
                0170      &       - 0.006758 _d 0 * COS ( 2. _d 0 * ALPHA )
                0171      &       + 0.000907 _d 0 * SIN ( 2. _d 0 * ALPHA )
                0172      &       - 0.002697 _d 0 * COS ( 3. _d 0 * ALPHA )
                0173      &       + 0.001480 _d 0 * SIN ( 3. _d 0 * ALPHA )
a303c567c2 Jean*0174
                0175 C note: alternative formulas include
                0176 C   1) formula from aim_surf_bc.F, neglecting eccentricity:
9675e9700b Jean*0177 C         ALPHA= 2. _d 0*PI*(TYEAR+10. _d 0/365. _d 0)
                0178 C         DECLI = COS(ALPHA) * ( -23.45 _d 0 * deg2rad)
a303c567c2 Jean*0179 C   2) formulas that accounts for minor astronomic effects, e.g.
                0180 C    Yallop, B. D., Position of the sun to 1 minute of arc precision,
                0181 C     H. M. Nautical Almanac Office, Royal Greenwich Observatory,
                0182 C     Herstmonceux Castle, Hailsham, Sussex BN27 1RP, 1977.
                0183
9675e9700b Jean*0184           ZC = COS(DECLI)
                0185           ZS = SIN(DECLI)
                0186           DO j = 1,sNy
                0187            DO i = 1,sNx
                0188
                0189             SJ = SIN(yC(i,j,bi,bj) * deg2rad)
                0190             CJ = COS(yC(i,j,bi,bj) * deg2rad)
                0191             TMPA = SJ*ZS
                0192             TMPB = CJ*ZC
94390f4f16 Gael*0193
9675e9700b Jean*0194             IF ( select_ZenAlbedo.EQ.3 ) THEN
a303c567c2 Jean*0195 C determine DAILY VARYING cos of solar zenith angle CZEN
                0196 C ------------------------------------------------------
                0197 C       (formula from Hartmann textbook, classic trigo)
9675e9700b Jean*0198              CZENdiurnal = TMPA + TMPB *
                0199      &         COS( 2. _d 0 *PI* TDAY + xC(i,j,bi,bj) * deg2rad )
a303c567c2 Jean*0200 C note: a more complicated hour angle formula is given by Yallop 1977
9675e9700b Jean*0201              IF ( CZENdiurnal .LE.0 ) CZENdiurnal = 0. _d 0
                0202              CZEN = CZENdiurnal
94390f4f16 Gael*0203
9675e9700b Jean*0204             ELSEIF ( select_ZenAlbedo.EQ.2 ) THEN
a303c567c2 Jean*0205 C determine DAILY MEAN cos of solar zenith angle CZEN
                0206 C ---------------------------------------------------
                0207 C       ( formula from aim_surf_bc.F <--> mean(CZEN*CZEN)/mean(CZEN) )
9675e9700b Jean*0208              TMPL = -TMPA/TMPB
                0209              IF (TMPL .GE. 1.0 _d 0) THEN
                0210               CZEN = 0.0 _d 0
                0211              ELSEIF (TMPL .LE. -1.0 _d 0) THEN
                0212               CZEN = (2.0 _d 0)*TMPA*PI
                0213               CZEN2= PI*((2.0 _d 0)*TMPA*TMPA + TMPB*TMPB)
                0214               CZEN = CZEN2/CZEN
                0215              ELSE
                0216               hlim = ACOS(TMPL)
                0217               CZEN = 2.0 _d 0*(TMPA*hlim + TMPB*SIN(hlim))
                0218               CZEN2= 2.0 _d 0*TMPA*TMPA*hlim
94390f4f16 Gael*0219      &          + 4.0 _d 0*TMPA*TMPB*SIN(hlim)
                0220      &          + TMPB*TMPB*( hlim + 0.5 _d 0*SIN(2.0 _d 0*hlim) )
9675e9700b Jean*0221               CZEN = CZEN2/CZEN
                0222              ENDIF
                0223              CZENdaily = CZEN
                0224 c            CZEN = CZENdaily
94390f4f16 Gael*0225
9675e9700b Jean*0226             ELSE
                0227               print *, 'select_ZenAlbedo is out of range'
                0228               STOP 'ABNORMAL END: S/R EXF_ZENITHANGLE'
                0229             ENDIF
94390f4f16 Gael*0230
a303c567c2 Jean*0231 C determine direct ocean albedo
                0232 C -----------------------------
                0233 C     (formula from Briegleb, Minnis, et al 1986)
                0234 C     comments on select_ZenAlbedo.GT.1 methods:
                0235 C     - CZENdaily as computed in aim was found to imply sizable biases in
                0236 C       daily upward sw fluxes.  It is not advised to use it, but it is kept
                0237 C       in connection to pkg/aim_v23.
                0238 C     - CZENdiurnal should never be used with daily mean input fields.
                0239 C       Furthermore, at this point, it is not advised to use it even with 6
                0240 C       hourly swdown input fields. This is because we simply time interpolate
                0241 C       between 6 hourly swdown fields, so each day there will be times when
                0242 C       CZENdiurnal correctly reflects that it is night time, but swdown.NE.0.
                0243 C       does not. CZENdiurnal may actually be rather harmful in this context,
                0244 C       since an inconsistency of phase between CZENdiurnal and swdown will
                0245 C       yield biases in daily mean upward sw fluxes. So ...
                0246
9675e9700b Jean*0247             ALBSEA1 =
                0248      &            ( ( 2.6 _d 0 / (CZEN**(1.7 _d 0) + 0.065 _d 0) )
94390f4f16 Gael*0249      &          + ( 15. _d 0 * (CZEN-0.1 _d 0) * (CZEN-0.5 _d 0)
                0250      &          * (CZEN-1.0 _d 0) ) ) / 100.0 _d 0
f71b2ceff6 Gael*0251
a303c567c2 Jean*0252 C determine overall albedo
                0253 C ------------------------
                0254 C       (approximation: half direct and half diffu.)
9675e9700b Jean*0255             zen_albedo (i,j,bi,bj) =
94390f4f16 Gael*0256      &          0.5 _d 0 * exf_albedo + 0.5 _d 0 * ALBSEA1
9675e9700b Jean*0257            ENDDO
d106b5e2d8 Gael*0258           ENDDO
9675e9700b Jean*0259
                0260 C        end if select_ZenAlbedo = 0, = 1, else
                0261          ENDIF
                0262
                0263 C       end bi,bj loops
d106b5e2d8 Gael*0264         ENDDO
                0265        ENDDO
                0266
a303c567c2 Jean*0267 C      end if ( useExfZenAlbedo )
9675e9700b Jean*0268       ENDIF
d106b5e2d8 Gael*0269
9675e9700b Jean*0270 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
d106b5e2d8 Gael*0271
                0272       IF ( useExfZenIncoming ) THEN
                0273
                0274        DO bj = myByLo(myThid),myByHi(myThid)
                0275         DO bi = myBxLo(myThid),myBxHi(myThid)
                0276
a303c567c2 Jean*0277 C compute incoming flux at the top of the atm.:
                0278 C ---------------------------------------------
                0279 C       (formula from Hartmann textbook, after Spencer 1971)
9675e9700b Jean*0280          ALPHA= 2. _d 0*PI*TYEAR
                0281          DECLI = 0.006918 _d 0
                0282      &         - 0.399912 _d 0 * COS ( 1. _d 0 * ALPHA )
                0283      &         + 0.070257 _d 0 * SIN ( 1. _d 0 * ALPHA )
                0284      &         - 0.006758 _d 0 * COS ( 2. _d 0 * ALPHA )
                0285      &         + 0.000907 _d 0 * SIN ( 2. _d 0 * ALPHA )
                0286      &         - 0.002697 _d 0 * COS ( 3. _d 0 * ALPHA )
                0287      &         + 0.001480 _d 0 * SIN ( 3. _d 0 * ALPHA )
                0288          dD0dDsq = 1.000110 _d 0
                0289      &         + 0.034221 _d 0 * COS ( 1. _d 0 * ALPHA )
                0290      &         + 0.001280 _d 0 * SIN ( 1. _d 0 * ALPHA )
                0291      &         + 0.000719 _d 0 * COS ( 2. _d 0 * ALPHA )
                0292      &         + 0.000077 _d 0 * SIN ( 2. _d 0 * ALPHA )
                0293          ZC = COS(DECLI)
                0294          ZS = SIN(DECLI)
                0295
                0296          DO j = 1,sNy
                0297           DO i = 1,sNx
                0298            SJ = SIN(yC(i,j,bi,bj) * deg2rad)
                0299            CJ = COS(yC(i,j,bi,bj) * deg2rad)
                0300            TMPA = SJ*ZS
                0301            TMPB = CJ*ZC
a303c567c2 Jean*0302 C DAILY VARYING value:
9675e9700b Jean*0303            CZEN = TMPA + TMPB *
                0304      &         COS( 2. _d 0 *PI*TDAY + xC(i,j,bi,bj)*deg2rad )
                0305            IF ( CZEN .LE.0 ) CZEN = 0. _d 0
                0306            FSOL = SOLC * dD0dDsq * MAX( 0. _d 0, CZEN )
                0307            zen_fsol_diurnal(i,j,bi,bj) = FSOL
                0308
a303c567c2 Jean*0309 C DAILY MEAN value:
9675e9700b Jean*0310            H0 = -TAN( yC(i,j,bi,bj) *deg2rad ) * TAN( DECLI )
                0311            IF ( H0.LT.-1. _d 0 ) H0 = -1. _d 0
                0312            IF ( H0.GT. 1. _d 0 ) H0 =  1. _d 0
                0313            H0 = ACOS( H0 )
                0314            FSOL = SOLC * dD0dDsq / PI
                0315      &          * ( H0 * TMPA + SIN(H0) * TMPB )
                0316            zen_fsol_daily(i,j,bi,bj) = FSOL
94390f4f16 Gael*0317
a303c567c2 Jean*0318 C note: an alternative for the DAILY MEAN is, as done in pkg/aim_v23,
9675e9700b Jean*0319 C          ALPHA= 2. _d 0*PI*(TYEAR+10. _d 0/365. _d 0)
                0320 C          CSR1=-0.796 _d 0*COS(ALPHA)
                0321 C          CSR2= 0.147 _d 0*COS(2. _d 0*ALPHA)-0.477 _d 0
                0322 C          FLAT2 = 1.5 _d 0*SJ**2 - 0.5 _d 0
                0323 C          FSOL = 0.25 _d 0 * SOLC
                0324 C    &          * MAX( 0. _d 0, 1. _d 0+CSR1*SJ+CSR2*FLAT2 )
                0325 C          zen_fsol_daily (i,j,bi,bj) = FSOL
94390f4f16 Gael*0326           ENDDO
                0327          ENDDO
9675e9700b Jean*0328
                0329 C       end bi,bj loops
94390f4f16 Gael*0330         ENDDO
                0331        ENDDO
                0332
a303c567c2 Jean*0333 C      end if ( useExfZenIncoming )
9675e9700b Jean*0334       ENDIF
d106b5e2d8 Gael*0335
94390f4f16 Gael*0336 #endif /* ALLOW_ZENITHANGLE */
                0337 #endif /* ALLOW_DOWNWARD_RADIATION */
                0338
                0339       RETURN
                0340       END