File indexing completed on 2018-03-02 18:40:02 UTC
view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
94390f4f16 Gael*0001 #include "EXF_OPTIONS.h"
0002
9675e9700b Jean*0003
0004
0005
0006 SUBROUTINE EXF_ZENITHANGLE(
0007 I myTime, myIter, myThid )
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
94390f4f16 Gael*0018 IMPLICIT NONE
0019
0020
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
94390f4f16 Gael*0033 _RL myTime
0034 INTEGER myIter
0035 INTEGER myThid
0036
0037 #ifdef ALLOW_DOWNWARD_RADIATION
0038 #ifdef ALLOW_ZENITHANGLE
9675e9700b Jean*0039
0040 #ifdef ALLOW_CAL
0041 INTEGER cal_IsLeap
0042 EXTERNAL cal_IsLeap
0043 #endif
0044
0045
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
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
94390f4f16 Gael*0061
a303c567c2 Jean*0062
0063
94390f4f16 Gael*0064 SOLC = 1368. _d 0
a303c567c2 Jean*0065
94390f4f16 Gael*0066
a303c567c2 Jean*0067
0068
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
0125
0126
0127
0128
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
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
0161 ELSE
94390f4f16 Gael*0162
a303c567c2 Jean*0163
0164
0165
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
0176
9675e9700b Jean*0177
0178
a303c567c2 Jean*0179
0180
0181
0182
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
0196
0197
9675e9700b Jean*0198 CZENdiurnal = TMPA + TMPB *
0199 & COS( 2. _d 0 *PI* TDAY + xC(i,j,bi,bj) * deg2rad )
a303c567c2 Jean*0200
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
0206
0207
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
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
0232
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244
0245
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
0253
0254
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
0261 ENDIF
0262
0263
d106b5e2d8 Gael*0264 ENDDO
0265 ENDDO
0266
a303c567c2 Jean*0267
9675e9700b Jean*0268 ENDIF
d106b5e2d8 Gael*0269
9675e9700b Jean*0270
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
0278
0279
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
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
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
9675e9700b Jean*0319
0320
0321
0322
0323
0324
0325
94390f4f16 Gael*0326 ENDDO
0327 ENDDO
9675e9700b Jean*0328
0329
94390f4f16 Gael*0330 ENDDO
0331 ENDDO
0332
a303c567c2 Jean*0333
9675e9700b Jean*0334 ENDIF
d106b5e2d8 Gael*0335
94390f4f16 Gael*0336 #endif /* ALLOW_ZENITHANGLE */
0337 #endif /* ALLOW_DOWNWARD_RADIATION */
0338
0339 RETURN
0340 END