| 1 |
C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_swrad.F,v 1.3 2004/06/17 21:45:29 molod Exp $ |
| 2 |
C $Name: $ |
| 3 |
|
| 4 |
#include "CPP_OPTIONS.h" |
| 5 |
subroutine swrio (nymd,nhms,bi,bj,ndswr,myid,istrip,npcs, |
| 6 |
. pz,tz,qz,pkht,oz,co2, |
| 7 |
. albvisdr,albvisdf,albnirdr,albnirdf, |
| 8 |
. dtradsw,dtswclr,radswg,swgclr, |
| 9 |
. fdifpar,fdirpar,osr,osrclr, |
| 10 |
. im,jm,lm,sige,sig,dsig,ptop, |
| 11 |
. nswcld,cldsw,cswmo,nswlz,swlz, |
| 12 |
. lpnt,imstturb,qliqave,fccave,landtype,xlats,xlons) |
| 13 |
|
| 14 |
implicit none |
| 15 |
#ifdef ALLOW_DIAGNOSTICS |
| 16 |
#include "diagnostics.h" |
| 17 |
#endif |
| 18 |
|
| 19 |
c Input Variables |
| 20 |
c --------------- |
| 21 |
integer nymd,nhms,ndswr,istrip,npcs,bi,bj |
| 22 |
|
| 23 |
integer im,jm,lm |
| 24 |
real ptop |
| 25 |
real sige(lm+1) |
| 26 |
real sig(lm) |
| 27 |
real dsig(lm) |
| 28 |
|
| 29 |
real pz(im,jm) |
| 30 |
real tz(im,jm,lm) |
| 31 |
real pkht(im,jm,lm) |
| 32 |
|
| 33 |
real co2 |
| 34 |
real oz(im,jm,lm) |
| 35 |
real qz(im,jm,lm) |
| 36 |
|
| 37 |
real albvisdr(im,jm) |
| 38 |
real albvisdf(im,jm) |
| 39 |
real albnirdr(im,jm) |
| 40 |
real albnirdf(im,jm) |
| 41 |
|
| 42 |
real radswg(im,jm) |
| 43 |
real swgclr(im,jm) |
| 44 |
real fdifpar(im,jm) |
| 45 |
real fdirpar(im,jm) |
| 46 |
real osr(im,jm) |
| 47 |
real osrclr(im,jm) |
| 48 |
real dtradsw(im,jm,lm) |
| 49 |
real dtswclr(im,jm,lm) |
| 50 |
|
| 51 |
integer nswcld,nswlz |
| 52 |
real cldsw(im,jm,lm) |
| 53 |
real cswmo(im,jm,lm) |
| 54 |
real swlz(im,jm,lm) |
| 55 |
|
| 56 |
logical lpnt |
| 57 |
integer imstturb |
| 58 |
real qliqave(im,jm,lm) |
| 59 |
real fccave(im,jm,lm) |
| 60 |
|
| 61 |
integer landtype(im,jm) |
| 62 |
|
| 63 |
c Local Variables |
| 64 |
c --------------- |
| 65 |
integer i,j,L,nn,nsecf,mid_level, low_level |
| 66 |
integer nb2,ntmstp,nymd2,nhms2 |
| 67 |
real getcon,grav,cp,undef,pcheck |
| 68 |
real ra,alf,reffw,reffi,tminv |
| 69 |
|
| 70 |
parameter ( reffw = 10.0 ) |
| 71 |
parameter ( reffi = 65.0 ) |
| 72 |
|
| 73 |
real alat(im,jm) |
| 74 |
real alon(im,jm) |
| 75 |
|
| 76 |
real PKZ(im,jm,lm) |
| 77 |
real PLZ(im,jm,lm) |
| 78 |
real tdry(im,jm,lm) |
| 79 |
real PLZE(im,jm,lm+1) |
| 80 |
real TEMP1(im,jm) |
| 81 |
real TEMP2(im,jm) |
| 82 |
real zenith (im,jm) |
| 83 |
real cldtot (im,jm,lm) |
| 84 |
real cldmxo (im,jm,lm) |
| 85 |
real totcld (im,jm) |
| 86 |
real cldlow (im,jm) |
| 87 |
real cldmid (im,jm) |
| 88 |
real cldhi (im,jm) |
| 89 |
real taulow (im,jm) |
| 90 |
real taumid (im,jm) |
| 91 |
real tauhi (im,jm) |
| 92 |
real tautype(im,jm,lm,3) |
| 93 |
real tau (im,jm,lm) |
| 94 |
real albedo(im,jm) |
| 95 |
|
| 96 |
real PK(ISTRIP,lm) |
| 97 |
real qzl(ISTRIP,lm), CLRO(ISTRIP,lm) |
| 98 |
real TZL(ISTRIP,lm) |
| 99 |
real OZL(ISTRIP,lm) |
| 100 |
real PLE(ISTRIP,lm+1) |
| 101 |
real COSZ(ISTRIP) |
| 102 |
|
| 103 |
real albuvdr(ISTRIP),albuvdf(ISTRIP) |
| 104 |
real albirdr(ISTRIP),albirdf(ISTRIP) |
| 105 |
real difpar (ISTRIP),dirpar (ISTRIP) |
| 106 |
|
| 107 |
real fdirir(istrip),fdifir(istrip) |
| 108 |
real fdiruv(istrip),fdifuv(istrip) |
| 109 |
|
| 110 |
real flux (istrip,lm+1) |
| 111 |
real fluxclr(istrip,lm+1) |
| 112 |
real dtsw (istrip,lm) |
| 113 |
real dtswc (istrip,lm) |
| 114 |
|
| 115 |
real taul (istrip,lm) |
| 116 |
real reff (istrip,lm,2) |
| 117 |
real tauc (istrip,lm,2) |
| 118 |
real taua (istrip,lm) |
| 119 |
real tstrip (istrip) |
| 120 |
|
| 121 |
logical first |
| 122 |
data first /.true./ |
| 123 |
|
| 124 |
integer koz, kh2o |
| 125 |
data KOZ /20/ |
| 126 |
data kh2o /18/ |
| 127 |
|
| 128 |
C ********************************************************************** |
| 129 |
C **** INITIALIZATION **** |
| 130 |
C ********************************************************************** |
| 131 |
|
| 132 |
grav = getcon('GRAVITY') |
| 133 |
cp = getcon('CP') |
| 134 |
undef = getcon('UNDEF') |
| 135 |
|
| 136 |
NTMSTP = nsecf(NDSWR) |
| 137 |
TMINV = 1./float(ntmstp) |
| 138 |
|
| 139 |
do j = 1,jm |
| 140 |
do i = 1,im |
| 141 |
PLZE(I,j,1) = SIGE(1)*PZ(I,j) + PTOP |
| 142 |
enddo |
| 143 |
enddo |
| 144 |
|
| 145 |
DO L = 1,lm |
| 146 |
do j = 1,jm |
| 147 |
DO I = 1,im |
| 148 |
PLZ (I,j,L ) = SIG (L) *PZ(I,j) + PTOP |
| 149 |
PLZE(I,j,L+1) = SIGE(L+1)*PZ(I,j) + PTOP |
| 150 |
ENDDO |
| 151 |
ENDDO |
| 152 |
ENDDO |
| 153 |
|
| 154 |
call pkappa ( pz,pkht,pkz,ptop,sige,dsig,im,jm,lm ) |
| 155 |
|
| 156 |
C Compute Temperature from Theta |
| 157 |
C ------------------------------ |
| 158 |
do L=1,lm |
| 159 |
do j=1,jm |
| 160 |
do i=1,im |
| 161 |
tdry(i,j,L) = tz(i,j,L)*pkz(i,j,L) |
| 162 |
enddo |
| 163 |
enddo |
| 164 |
enddo |
| 165 |
|
| 166 |
c Determine Level Indices for Low-Mid-High Cloud Regions |
| 167 |
c ------------------------------------------------------ |
| 168 |
low_level = lm |
| 169 |
mid_level = lm |
| 170 |
do L = lm-1,1,-1 |
| 171 |
pcheck = (1000.-ptop)*sig(l) + ptop |
| 172 |
if (pcheck.gt.700.0) low_level = L |
| 173 |
if (pcheck.gt.400.0) mid_level = L |
| 174 |
enddo |
| 175 |
|
| 176 |
if (first .and. myid.eq.0 ) then |
| 177 |
print * |
| 178 |
print *,'Low-Level Clouds are Grouped between levels: ', |
| 179 |
. lm,' and ',low_level |
| 180 |
print *,'Mid-Level Clouds are Grouped between levels: ', |
| 181 |
. low_level-1,' and ',mid_level |
| 182 |
print * |
| 183 |
first = .false. |
| 184 |
endif |
| 185 |
|
| 186 |
C ********************************************************************** |
| 187 |
C **** CALCULATE COSINE OF THE ZENITH ANGLE **** |
| 188 |
C ********************************************************************** |
| 189 |
|
| 190 |
CALL ASTRO ( NYMD, NHMS, ALAT,ALON, im*jm, TEMP1,RA ) |
| 191 |
NYMD2 = NYMD |
| 192 |
NHMS2 = NHMS |
| 193 |
CALL TICK ( NYMD2, NHMS2, NTMSTP ) |
| 194 |
CALL ASTRO ( NYMD2, NHMS2, ALAT,ALON, im*jm, TEMP2,RA ) |
| 195 |
|
| 196 |
do j = 1,jm |
| 197 |
do i = 1,im |
| 198 |
zenith(I,j) = TEMP1(I,j) + TEMP2(I,j) |
| 199 |
IF( zenith(I,j) .GT. 0.0 ) THEN |
| 200 |
zenith(I,j) = (2./3.)*( zenith(I,j)-TEMP2(I,j)*TEMP1(I,j) |
| 201 |
. / zenith(I,j) ) |
| 202 |
ENDIF |
| 203 |
ENDDO |
| 204 |
ENDDO |
| 205 |
|
| 206 |
|
| 207 |
C ********************************************************************** |
| 208 |
c **** Compute Two-Dimension Total Cloud Fraction (0-1) **** |
| 209 |
C ********************************************************************** |
| 210 |
|
| 211 |
c Initialize Clear Sky Probability for Low, Mid, and High Clouds |
| 212 |
c -------------------------------------------------------------- |
| 213 |
do j =1,jm |
| 214 |
do i =1,im |
| 215 |
cldlow(i,j) = 0.0 |
| 216 |
cldmid(i,j) = 0.0 |
| 217 |
cldhi (i,j) = 0.0 |
| 218 |
enddo |
| 219 |
enddo |
| 220 |
|
| 221 |
c Adjust cloud fractions and cloud liquid water due to moist turbulence |
| 222 |
c --------------------------------------------------------------------- |
| 223 |
if(imstturb.ne.0) then |
| 224 |
do L =1,lm |
| 225 |
do j =1,jm |
| 226 |
do i =1,im |
| 227 |
cldtot(i,j,L)=min(1.0,max(cldsw(i,j,L),fccave(i,j,L)/imstturb)) |
| 228 |
cldmxo(i,j,L)=min(1.0,cswmo(i,j,L)) |
| 229 |
swlz(i,j,L)=swlz(i,j,L)+qliqave(i,j,L)/imstturb |
| 230 |
enddo |
| 231 |
enddo |
| 232 |
enddo |
| 233 |
else |
| 234 |
do L =1,lm |
| 235 |
do j =1,jm |
| 236 |
do i =1,im |
| 237 |
cldtot(i,j,L) = min( 1.0,cldsw(i,j,L) ) |
| 238 |
cldmxo(i,j,L) = min( 1.0,cswmo(i,j,L) ) |
| 239 |
enddo |
| 240 |
enddo |
| 241 |
enddo |
| 242 |
endif |
| 243 |
|
| 244 |
c Compute 3-D Cloud Fractions |
| 245 |
c --------------------------- |
| 246 |
if( nswcld.ne.0 ) then |
| 247 |
do L = 1,lm |
| 248 |
do j = 1,jm |
| 249 |
do i = 1,im |
| 250 |
c Compute Low-Mid-High Maximum Overlap Cloud Fractions |
| 251 |
c ---------------------------------------------------- |
| 252 |
if( L.lt.mid_level ) then |
| 253 |
cldhi (i,j) = max( cldtot(i,j,L),cldhi (i,j) ) |
| 254 |
elseif( L.ge.mid_level .and. |
| 255 |
. L.lt.low_level ) then |
| 256 |
cldmid(i,j) = max( cldtot(i,j,L),cldmid(i,j) ) |
| 257 |
elseif( L.ge.low_level ) then |
| 258 |
cldlow(i,j) = max( cldtot(i,j,L),cldlow(i,j) ) |
| 259 |
endif |
| 260 |
|
| 261 |
enddo |
| 262 |
enddo |
| 263 |
enddo |
| 264 |
endif |
| 265 |
|
| 266 |
c Totcld => Product of Clear Sky Probabilities for Low, Mid, and High Clouds |
| 267 |
c -------------------------------------------------------------------------- |
| 268 |
do j = 1,jm |
| 269 |
do i = 1,im |
| 270 |
totcld(i,j) = 1.0 - (1.-cldhi (i,j)) |
| 271 |
. * (1.-cldmid(i,j)) |
| 272 |
. * (1.-cldlow(i,j)) |
| 273 |
enddo |
| 274 |
enddo |
| 275 |
|
| 276 |
c Compute Cloud Diagnostics |
| 277 |
c ------------------------- |
| 278 |
if(icldfrc.gt.0) then |
| 279 |
do j=1,jm |
| 280 |
do i=1,im |
| 281 |
qdiag(i,j,icldfrc,bi,bj) = qdiag(i,j,icldfrc,bi,bj) + totcld(i,j) |
| 282 |
enddo |
| 283 |
enddo |
| 284 |
ncldfrc = ncldfrc + 1 |
| 285 |
endif |
| 286 |
|
| 287 |
if( icldras.gt.0 ) then |
| 288 |
do L=1,lm |
| 289 |
do j=1,jm |
| 290 |
do i=1,im |
| 291 |
qdiag(i,j,icldras+L-1,bi,bj) = qdiag(i,j,icldras+L-1,bi,bj) + |
| 292 |
. cswmo(i,j,L) |
| 293 |
enddo |
| 294 |
enddo |
| 295 |
enddo |
| 296 |
ncldras = ncldras + 1 |
| 297 |
endif |
| 298 |
|
| 299 |
if( icldtot.gt.0 ) then |
| 300 |
do L=1,lm |
| 301 |
do j=1,jm |
| 302 |
do i=1,im |
| 303 |
qdiag(i,j,icldtot+L-1,bi,bj) = qdiag(i,j,icldtot+L-1,bi,bj) + |
| 304 |
. cldtot(i,j,L) |
| 305 |
enddo |
| 306 |
enddo |
| 307 |
enddo |
| 308 |
ncldtot = ncldtot + 1 |
| 309 |
endif |
| 310 |
|
| 311 |
if( icldlow.gt.0 ) then |
| 312 |
do j=1,jm |
| 313 |
do i=1,im |
| 314 |
qdiag(i,j,icldlow,bi,bj) = qdiag(i,j,icldlow,bi,bj) + cldlow(i,j) |
| 315 |
enddo |
| 316 |
enddo |
| 317 |
ncldlow = ncldlow + 1 |
| 318 |
endif |
| 319 |
|
| 320 |
if( icldmid.gt.0 ) then |
| 321 |
do j=1,jm |
| 322 |
do i=1,im |
| 323 |
qdiag(i,j,icldmid,bi,bj) = qdiag(i,j,icldmid,bi,bj) + cldmid(i,j) |
| 324 |
enddo |
| 325 |
enddo |
| 326 |
ncldmid = ncldmid + 1 |
| 327 |
endif |
| 328 |
|
| 329 |
if( icldhi.gt.0 ) then |
| 330 |
do j=1,jm |
| 331 |
do i=1,im |
| 332 |
qdiag(i,j,icldhi,bi,bj) = qdiag(i,j,icldhi,bi,bj) + cldhi(i,j) |
| 333 |
enddo |
| 334 |
enddo |
| 335 |
ncldhi = ncldhi + 1 |
| 336 |
endif |
| 337 |
|
| 338 |
if( ilzrad.gt.0 ) then |
| 339 |
do L=1,lm |
| 340 |
do j=1,jm |
| 341 |
do i=1,im |
| 342 |
qdiag(i,j,ilzrad+L-1,bi,bj) = qdiag(i,j,ilzrad+L-1,bi,bj) + |
| 343 |
. swlz(i,j,L)*1.0e6 |
| 344 |
enddo |
| 345 |
enddo |
| 346 |
enddo |
| 347 |
nlzrad = nlzrad + 1 |
| 348 |
endif |
| 349 |
|
| 350 |
c Albedo Diagnostics |
| 351 |
c ------------------ |
| 352 |
if( ialbvisdr.gt.0 ) then |
| 353 |
do j=1,jm |
| 354 |
do i=1,im |
| 355 |
qdiag(i,j,ialbvisdr,bi,bj) = qdiag(i,j,ialbvisdr,bi,bj) + |
| 356 |
. albvisdr(i,j) |
| 357 |
enddo |
| 358 |
enddo |
| 359 |
nalbvisdr = nalbvisdr + 1 |
| 360 |
endif |
| 361 |
|
| 362 |
if( ialbvisdf.gt.0 ) then |
| 363 |
do j=1,jm |
| 364 |
do i=1,im |
| 365 |
qdiag(i,j,ialbvisdf,bi,bj) = qdiag(i,j,ialbvisdf,bi,bj) + |
| 366 |
. albvisdf(i,j) |
| 367 |
enddo |
| 368 |
enddo |
| 369 |
nalbvisdf = nalbvisdf + 1 |
| 370 |
endif |
| 371 |
|
| 372 |
if( ialbnirdr.gt.0 ) then |
| 373 |
do j=1,jm |
| 374 |
do i=1,im |
| 375 |
qdiag(i,j,ialbnirdr,bi,bj) = qdiag(i,j,ialbnirdr,bi,bj) + |
| 376 |
. albnirdr(i,j) |
| 377 |
enddo |
| 378 |
enddo |
| 379 |
nalbnirdr = nalbnirdr + 1 |
| 380 |
endif |
| 381 |
|
| 382 |
if( ialbnirdf.gt.0 ) then |
| 383 |
do j=1,jm |
| 384 |
do i=1,im |
| 385 |
qdiag(i,j,ialbnirdf,bi,bj) = qdiag(i,j,ialbnirdf,bi,bj) + |
| 386 |
. albnirdf(i,j) |
| 387 |
enddo |
| 388 |
enddo |
| 389 |
nalbnirdf = nalbnirdf + 1 |
| 390 |
endif |
| 391 |
|
| 392 |
C Compute Optical Thicknesses and Diagnostics |
| 393 |
C ------------------------------------------- |
| 394 |
call opthk(tdry,plz,plze,swlz,cldtot,cldmxo,landtype,im,jm,lm, |
| 395 |
. tautype) |
| 396 |
|
| 397 |
do L = 1,lm |
| 398 |
do j = 1,jm |
| 399 |
do i = 1,im |
| 400 |
tau(i,j,L) = tautype(i,j,L,1)+tautype(i,j,L,2)+tautype(i,j,L,3) |
| 401 |
enddo |
| 402 |
enddo |
| 403 |
enddo |
| 404 |
|
| 405 |
if( itauave.gt.0 ) then |
| 406 |
do L=1,lm |
| 407 |
do j=1,jm |
| 408 |
do i=1,im |
| 409 |
qdiag(i,j,itauave+L-1,bi,bj) = qdiag(i,j,itauave+L-1,bi,bj) + |
| 410 |
. tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L)) |
| 411 |
enddo |
| 412 |
enddo |
| 413 |
enddo |
| 414 |
ntauave = ntauave + 1 |
| 415 |
endif |
| 416 |
|
| 417 |
if( itaucld.gt.0 ) then |
| 418 |
do L=1,lm |
| 419 |
do j=1,jm |
| 420 |
do i=1,im |
| 421 |
if( cldtot(i,j,L).ne.0.0 ) then |
| 422 |
qdiag(i,j,itaucld +L-1,bi,bj) = qdiag(i,j,itaucld +L-1,bi,bj) + |
| 423 |
. tau(i,j,L)*100/(plze(i,j,L+1)-plze(i,j,L)) |
| 424 |
qdiag(i,j,itaucldc+L-1,bi,bj) = |
| 425 |
. qdiag(i,j,itaucldc+L-1,bi,bj) + 1.0 |
| 426 |
endif |
| 427 |
enddo |
| 428 |
enddo |
| 429 |
enddo |
| 430 |
endif |
| 431 |
|
| 432 |
c Compute Low, Mid, and High Cloud Optical Depth Diagnostics |
| 433 |
c ---------------------------------------------------------- |
| 434 |
if( itaulow.ne.0 ) then |
| 435 |
do j = 1,jm |
| 436 |
do i = 1,im |
| 437 |
if( cldlow(i,j).ne.0.0 ) then |
| 438 |
taulow(i,j) = 0.0 |
| 439 |
do L = low_level,lm |
| 440 |
taulow(i,j) = taulow(i,j) + tau(i,j,L) |
| 441 |
enddo |
| 442 |
qdiag(i,j,itaulow,bi,bj ) = qdiag(i,j,itaulow,bi,bj ) + |
| 443 |
. taulow(i,j) |
| 444 |
qdiag(i,j,itaulowc,bi,bj) = qdiag(i,j,itaulowc,bi,bj) + 1.0 |
| 445 |
endif |
| 446 |
enddo |
| 447 |
enddo |
| 448 |
endif |
| 449 |
|
| 450 |
if( itaumid.ne.0 ) then |
| 451 |
do j = 1,jm |
| 452 |
do i = 1,im |
| 453 |
if( cldmid(i,j).ne.0.0 ) then |
| 454 |
taumid(i,j) = 0.0 |
| 455 |
do L = mid_level,low_level+1 |
| 456 |
taumid(i,j) = taumid(i,j) + tau(i,j,L) |
| 457 |
enddo |
| 458 |
qdiag(i,j,itaumid,bi,bj ) = qdiag(i,j,itaumid,bi,bj ) + |
| 459 |
. taumid(i,j) |
| 460 |
qdiag(i,j,itaumidc,bi,bj) = qdiag(i,j,itaumidc,bi,bj) + 1.0 |
| 461 |
endif |
| 462 |
enddo |
| 463 |
enddo |
| 464 |
endif |
| 465 |
|
| 466 |
if( itauhi.ne.0 ) then |
| 467 |
do j = 1,jm |
| 468 |
do i = 1,im |
| 469 |
if( cldhi(i,j).ne.0.0 ) then |
| 470 |
tauhi(i,j) = 0.0 |
| 471 |
do L = 1,mid_level+1 |
| 472 |
tauhi(i,j) = tauhi(i,j) + tau(i,j,L) |
| 473 |
enddo |
| 474 |
qdiag(i,j,itauhi,bi,bj ) = qdiag(i,j,itauhi,bi,bj ) + |
| 475 |
. tauhi(i,j) |
| 476 |
qdiag(i,j,itauhic,bi,bj) = qdiag(i,j,itauhic,bi,bj) + 1.0 |
| 477 |
endif |
| 478 |
enddo |
| 479 |
enddo |
| 480 |
endif |
| 481 |
|
| 482 |
C*********************************************************************** |
| 483 |
C **** LOOP OVER GLOBAL REGIONS **** |
| 484 |
C ********************************************************************** |
| 485 |
|
| 486 |
do 1000 nn = 1,npcs |
| 487 |
|
| 488 |
C ********************************************************************** |
| 489 |
C **** VARIABLE INITIALIZATION **** |
| 490 |
C ********************************************************************** |
| 491 |
|
| 492 |
CALL STRIP ( zenith,COSZ,im*jm,ISTRIP,1,NN ) |
| 493 |
|
| 494 |
CALL STRIP ( plze, ple ,im*jm,ISTRIP,lm+1,NN) |
| 495 |
CALL STRIP ( pkz , pk ,im*jm,ISTRIP,lm ,NN) |
| 496 |
CALL STRIP ( tdry, tzl ,im*jm,ISTRIP,lm ,NN) |
| 497 |
CALL STRIP ( qz , qzl ,im*jm,ISTRIP,lm ,NN) |
| 498 |
CALL STRIP ( oz , ozl ,im*jm,ISTRIP,lm ,NN) |
| 499 |
CALL STRIP ( tau , taul ,im*jm,ISTRIP,lm ,NN) |
| 500 |
|
| 501 |
CALL STRIP ( albvisdr,albuvdr,im*jm,ISTRIP,1,NN ) |
| 502 |
CALL STRIP ( albvisdf,albuvdf,im*jm,ISTRIP,1,NN ) |
| 503 |
CALL STRIP ( albnirdr,albirdr,im*jm,ISTRIP,1,NN ) |
| 504 |
CALL STRIP ( albnirdf,albirdf,im*jm,ISTRIP,1,NN ) |
| 505 |
|
| 506 |
call strip ( cldtot,clro,im*jm,istrip,lm,nn ) |
| 507 |
|
| 508 |
c Partition Tau between Water and Ice Particles |
| 509 |
c --------------------------------------------- |
| 510 |
do L= 1,lm |
| 511 |
do i= 1,istrip |
| 512 |
alf = min( max((tzl(i,l)-253.15)/20.,0.) ,1.) |
| 513 |
taua(i,L) = 0. |
| 514 |
|
| 515 |
if( alf.ne.0.0 .and. alf.ne.1.0 ) then |
| 516 |
tauc(i,L,1) = taul(i,L)/(1.+alf/(1-alf) * (reffi/reffw*0.80) ) |
| 517 |
tauc(i,L,2) = taul(i,L)/(1.+(1-alf)/alf * (reffw/reffi*1.25) ) |
| 518 |
endif |
| 519 |
|
| 520 |
if( alf.eq.0.0 ) then |
| 521 |
tauc(i,L,1) = taul(i,L) |
| 522 |
tauc(i,L,2) = 0.0 |
| 523 |
endif |
| 524 |
|
| 525 |
if( alf.eq.1.0 ) then |
| 526 |
tauc(i,L,1) = 0.0 |
| 527 |
tauc(i,L,2) = taul(i,L) |
| 528 |
endif |
| 529 |
|
| 530 |
reff(i,L,1) = reffi |
| 531 |
reff(i,L,2) = reffw |
| 532 |
enddo |
| 533 |
enddo |
| 534 |
|
| 535 |
call sorad ( istrip,1,1,lm,ple,tzl,qzl,ozl,co2, |
| 536 |
. tauc,reff,clro,mid_level,low_level,taua, |
| 537 |
. albirdr,albirdf,albuvdr,albuvdf,cosz, |
| 538 |
. flux,fluxclr,fdirir,fdifir,dirpar,difpar, |
| 539 |
. fdiruv,fdifuv ) |
| 540 |
|
| 541 |
C ********************************************************************** |
| 542 |
C **** Compute Mass-Weighted Theta Tendencies from Fluxes **** |
| 543 |
C ********************************************************************** |
| 544 |
|
| 545 |
do l=1,lm |
| 546 |
alf = grav/(cp*dsig(L)*100) |
| 547 |
do i=1,istrip |
| 548 |
dtsw (i,L) = alf*( flux (i,L)-flux (i,L+1) )/pk(i,L) |
| 549 |
dtswc(i,L) = alf*( fluxclr(i,L)-fluxclr(i,L+1) )/pk(i,L) |
| 550 |
enddo |
| 551 |
enddo |
| 552 |
|
| 553 |
call paste ( dtsw , dtradsw ,istrip,im*jm,lm,nn ) |
| 554 |
call paste ( dtswc, dtswclr ,istrip,im*jm,lm,nn ) |
| 555 |
|
| 556 |
call paste ( flux (1,1),osr ,istrip,im*jm,1,nn ) |
| 557 |
call paste ( fluxclr(1,1),osrclr,istrip,im*jm,1,nn ) |
| 558 |
|
| 559 |
call paste ( flux (1,lm+1),radswg,istrip,im*jm,1,nn ) |
| 560 |
call paste ( fluxclr(1,lm+1),swgclr,istrip,im*jm,1,nn ) |
| 561 |
|
| 562 |
call paste ( dirpar,fdirpar,istrip,im*jm,1,nn ) |
| 563 |
call paste ( difpar,fdifpar,istrip,im*jm,1,nn ) |
| 564 |
|
| 565 |
c Calculate Mean Albedo |
| 566 |
c --------------------- |
| 567 |
do i=1,istrip |
| 568 |
if( cosz(i).gt.0.0 ) then |
| 569 |
tstrip(i) = 1.0 - flux(i,lm+1)/ |
| 570 |
. ( fdirir(i)+fdifir(i)+dirpar(i)+difpar(i) + fdiruv(i)+fdifuv(i) ) |
| 571 |
if( tstrip(i).lt.0.0 ) tstrip(i) = undef |
| 572 |
else |
| 573 |
tstrip(i) = undef |
| 574 |
endif |
| 575 |
enddo |
| 576 |
call paste ( tstrip,albedo,istrip,im*jm,1,nn ) |
| 577 |
|
| 578 |
1000 CONTINUE |
| 579 |
|
| 580 |
c Mean Albedo Diagnostic |
| 581 |
c ---------------------- |
| 582 |
if( ialbedo.gt.0 ) then |
| 583 |
do j=1,jm |
| 584 |
do i=1,im |
| 585 |
if( albedo(i,j).ne.undef ) then |
| 586 |
qdiag(i,j,ialbedo,bi,bj ) = qdiag(i,j,ialbedo,bi,bj )+albedo(i,j) |
| 587 |
qdiag(i,j,ialbedoc,bi,bj) = qdiag(i,j,ialbedoc,bi,bj) + 1.0 |
| 588 |
endif |
| 589 |
enddo |
| 590 |
enddo |
| 591 |
endif |
| 592 |
|
| 593 |
C ********************************************************************** |
| 594 |
C **** ZERO-OUT OR BUMP COUNTERS **** |
| 595 |
C ********************************************************************** |
| 596 |
|
| 597 |
nswlz = 0 |
| 598 |
nswcld = 0 |
| 599 |
imstturb = 0 |
| 600 |
|
| 601 |
do L = 1,lm |
| 602 |
do j = 1,jm |
| 603 |
do i = 1,im |
| 604 |
fccave(i,j,L) = 0. |
| 605 |
qliqave(i,j,L) = 0. |
| 606 |
enddo |
| 607 |
enddo |
| 608 |
enddo |
| 609 |
|
| 610 |
return |
| 611 |
end |
| 612 |
subroutine opthk ( tl,pl,ple,lz,cf,cfm,lwi,im,jm,lm,tau ) |
| 613 |
C*********************************************************************** |
| 614 |
C |
| 615 |
C PURPOSE: |
| 616 |
C ======== |
| 617 |
C Compute cloud optical thickness using temperature and |
| 618 |
C cloud fractions. |
| 619 |
C |
| 620 |
C INPUT: |
| 621 |
C ====== |
| 622 |
C tl ......... Temperature at Model Layers (K) |
| 623 |
C pl ......... Model Layer Pressures (mb) |
| 624 |
C ple ........ Model Edge Pressures (mb) |
| 625 |
C lz ......... Diagnosed Convective and Large-Scale Cloud Water (g/g) |
| 626 |
C cf ......... Total Cloud Fraction (Random + Maximum Overlap) |
| 627 |
C cfm ........ Maximum Overlap Cloud Fraction |
| 628 |
C lwi ........ Surface Land Type |
| 629 |
C im ......... Longitudinal dimension |
| 630 |
C jm ......... Latitudinal dimension |
| 631 |
C lm ......... Vertical dimension |
| 632 |
C |
| 633 |
C OUTPUT: |
| 634 |
C ======= |
| 635 |
C tau ........ Cloud Optical Thickness (non-dimensional) |
| 636 |
C tau(im,jm,lm,1): Suspended Ice |
| 637 |
C tau(im,jm,lm,2): Suspended Water |
| 638 |
C tau(im,jm,lm,3): Raindrops |
| 639 |
C |
| 640 |
C*********************************************************************** |
| 641 |
C* GODDARD LABORATORY FOR ATMOSPHERES * |
| 642 |
C*********************************************************************** |
| 643 |
|
| 644 |
implicit none |
| 645 |
|
| 646 |
integer im,jm,lm,i,j,L |
| 647 |
|
| 648 |
real tl(im,jm,lm) |
| 649 |
real pl(im,jm,lm) |
| 650 |
real ple(im,jm,lm+1) |
| 651 |
real lz(im,jm,lm) |
| 652 |
real cf(im,jm,lm) |
| 653 |
real cfm(im,jm,lm) |
| 654 |
real tau(im,jm,lm,3) |
| 655 |
integer lwi(im,jm) |
| 656 |
|
| 657 |
real dp, alf, fracls, fraccu |
| 658 |
real tauice, tauh2o, tauras |
| 659 |
|
| 660 |
c Compute Cloud Optical Depths |
| 661 |
c ---------------------------- |
| 662 |
do L=1,lm |
| 663 |
do j=1,jm |
| 664 |
do i=1,im |
| 665 |
alf = min( max(( tl(i,j,L)-233.15)/20.,0.) ,1.) |
| 666 |
dp = ple(i,j,L+1)-ple(i,j,L) |
| 667 |
tau(i,j,L,1) = 0.0 |
| 668 |
tau(i,j,L,2) = 0.0 |
| 669 |
tau(i,j,L,3) = 0.0 |
| 670 |
if( cf(i,j,L).ne.0.0 ) then |
| 671 |
|
| 672 |
c Determine fraction of large-scale and cumulus clouds |
| 673 |
c ---------------------------------------------------- |
| 674 |
fracls = ( cf(i,j,L)-cfm(i,j,L) )/cf(i,j,L) |
| 675 |
fraccu = 1.0-fracls |
| 676 |
|
| 677 |
c Define tau for large-scale ice and water clouds |
| 678 |
c Note: tauice is scaled between (0.02 & .2) for: 0 < lz < 2 mg/kg over Land |
| 679 |
c Note: tauh2o is scaled between (0.20 & 20) for: 0 < lz < 5 mg/kg over Land |
| 680 |
c Note: tauh2o is scaled between (0.20 & 12) for: 0 < lz < 50 mg/kg over Ocean |
| 681 |
c ------------------------------------------------------------------------------- |
| 682 |
|
| 683 |
c Large-Scale Ice |
| 684 |
c --------------- |
| 685 |
tauice = max( 0.0002, 0.002*min(500*lz(i,j,L)*1000,1.0) ) |
| 686 |
tau(i,j,L,1) = fracls*(1-alf)*tauice*dp |
| 687 |
|
| 688 |
c Large-Scale Water |
| 689 |
c ----------------- |
| 690 |
C Over Land |
| 691 |
if( lwi(i,j).le.10 ) then |
| 692 |
tauh2o = max( 0.0020, 0.200*min(200*lz(i,j,L)*1000,1.0) ) |
| 693 |
tau(i,j,L,3) = fracls*alf*tauh2o*dp |
| 694 |
else |
| 695 |
C Non-Precipitation Clouds Over Ocean |
| 696 |
if( lz(i,j,L).eq.0.0 ) then |
| 697 |
tauh2o = .12 |
| 698 |
tau(i,j,L,2) = fracls*alf*tauh2o*dp |
| 699 |
else |
| 700 |
C Over Ocean |
| 701 |
tauh2o = max( 0.0020, 0.120*min( 20*lz(i,j,L)*1000,1.0) ) |
| 702 |
tau(i,j,L,3) = fracls*alf*tauh2o*dp |
| 703 |
endif |
| 704 |
endif |
| 705 |
|
| 706 |
c Sub-Grid Convective |
| 707 |
c ------------------- |
| 708 |
if( tl(i,j,L).gt.210.0 ) then |
| 709 |
tauras = .16 |
| 710 |
tau(i,j,L,3) = tau(i,j,L,3) + fraccu*tauras*dp |
| 711 |
else |
| 712 |
tauras = tauice |
| 713 |
tau(i,j,L,1) = tau(i,j,L,1) + fraccu*tauras*dp |
| 714 |
endif |
| 715 |
|
| 716 |
c Define tau for large-scale and cumulus clouds |
| 717 |
c --------------------------------------------- |
| 718 |
ccc tau(i,j,L) = ( fracls*( alf*tauh2o + (1-alf)*tauice ) |
| 719 |
ccc . + fraccu*tauras )*dp |
| 720 |
|
| 721 |
endif |
| 722 |
|
| 723 |
enddo |
| 724 |
enddo |
| 725 |
enddo |
| 726 |
|
| 727 |
return |
| 728 |
end |
| 729 |
subroutine sorad(m,n,ndim,np,pl,ta,wa,oa,co2, |
| 730 |
* taucld,reff,fcld,ict,icb, |
| 731 |
* taual,rsirbm,rsirdf,rsuvbm,rsuvdf,cosz, |
| 732 |
* flx,flc,fdirir,fdifir,fdirpar,fdifpar, |
| 733 |
* fdiruv,fdifuv) |
| 734 |
c******************************************************************** |
| 735 |
c |
| 736 |
c This routine computes solar fluxes due to the absoption by water |
| 737 |
c vapor, ozone, co2, o2, clouds, and aerosols and due to the |
| 738 |
c scattering by clouds, aerosols, and gases. |
| 739 |
c |
| 740 |
c This is a vectorized code. It computes the fluxes simultaneous for |
| 741 |
c (m x n) soundings, which is a subset of the (m x ndim) soundings. |
| 742 |
c In a global climate model, m and ndim correspond to the numbers of |
| 743 |
c grid boxes in the zonal and meridional directions, respectively. |
| 744 |
c |
| 745 |
c Ice and liquid cloud particles are allowed to co-exist in any of the |
| 746 |
c np layers. Two sets of cloud parameters are required as inputs, one |
| 747 |
c for ice paticles and the other for liquid particles. These parameters |
| 748 |
c are optical thickness (taucld) and effective particle size (reff). |
| 749 |
c |
| 750 |
c If no information is available for reff, a default value of |
| 751 |
c 10 micron for liquid water and 75 micron for ice can be used. |
| 752 |
c |
| 753 |
c Clouds are grouped into high, middle, and low clouds separated by the |
| 754 |
c level indices ict and icb. For detail, see the subroutine cldscale. |
| 755 |
c |
| 756 |
c----- Input parameters: |
| 757 |
c units size |
| 758 |
c number of soundings in zonal direction (m) n/d 1 |
| 759 |
c number of soundings in meridional direction (n) n/d 1 |
| 760 |
c maximum number of soundings in n/d 1 |
| 761 |
c meridional direction (ndim) |
| 762 |
c number of atmospheric layers (np) n/d 1 |
| 763 |
c level pressure (pl) mb m*ndim*(np+1) |
| 764 |
c layer temperature (ta) k m*ndim*np |
| 765 |
c layer specific humidity (wa) gm/gm m*ndim*np |
| 766 |
c layer ozone concentration (oa) gm/gm m*ndim*np |
| 767 |
c co2 mixing ratio by volumn (co2) parts/part 1 |
| 768 |
c cloud optical thickness (taucld) n/d m*ndim*np*2 |
| 769 |
c index 1 for ice particles |
| 770 |
c index 2 for liquid drops |
| 771 |
c effective cloud-particle size (reff) micrometer m*ndim*np*2 |
| 772 |
c index 1 for ice particles |
| 773 |
c index 2 for liquid drops |
| 774 |
c cloud amount (fcld) fraction m*ndim*np |
| 775 |
c level index separating high and middle n/d 1 |
| 776 |
c clouds (ict) |
| 777 |
c level index separating middle and low clouds n/d 1 |
| 778 |
c clouds (icb) |
| 779 |
c aerosol optical thickness (taual) n/d m*ndim*np |
| 780 |
c solar ir surface albedo for beam fraction m*ndim |
| 781 |
c radiation (rsirbm) |
| 782 |
c solar ir surface albedo for diffuse fraction m*ndim |
| 783 |
c radiation (rsirdf) |
| 784 |
c uv + par surface albedo for beam fraction m*ndim |
| 785 |
c radiation (rsuvbm) |
| 786 |
c uv + par surface albedo for diffuse fraction m*ndim |
| 787 |
c radiation (rsuvdf) |
| 788 |
c cosine of solar zenith angle (cosz) n/d m*ndim |
| 789 |
c |
| 790 |
c----- Output parameters |
| 791 |
c |
| 792 |
c all-sky flux (downward minus upward) (flx) fraction m*ndim*(np+1) |
| 793 |
c clear-sky flux (downward minus upward) (flc) fraction m*ndim*(np+1) |
| 794 |
c all-sky direct downward ir (0.7-10 micron) |
| 795 |
c flux at the surface (fdirir) fraction m*ndim |
| 796 |
c all-sky diffuse downward ir flux at |
| 797 |
c the surface (fdifir) fraction m*ndim |
| 798 |
c all-sky direct downward par (0.4-0.7 micron) |
| 799 |
c flux at the surface (fdirpar) fraction m*ndim |
| 800 |
c all-sky diffuse downward par flux at |
| 801 |
c the surface (fdifpar) fraction m*ndim |
| 802 |
* |
| 803 |
c----- Notes: |
| 804 |
c |
| 805 |
c (1) The unit of flux is fraction of the incoming solar radiation |
| 806 |
c at the top of the atmosphere. Therefore, fluxes should |
| 807 |
c be equal to flux multiplied by the extra-terrestrial solar |
| 808 |
c flux and the cosine of solar zenith angle. |
| 809 |
c (2) Clouds and aerosols can be included in any layers by specifying |
| 810 |
c fcld(i,j,k), taucld(i,j,k,*) and taual(i,j,k), k=1,np. |
| 811 |
c For an atmosphere without clouds and aerosols, |
| 812 |
c set fcld(i,j,k)=taucld(i,j,k,*)=taual(i,j,k)=0.0. |
| 813 |
c (3) Aerosol single scattering albedos and asymmetry |
| 814 |
c factors are specified in the subroutines solir and soluv. |
| 815 |
c (4) pl(i,j,1) is the pressure at the top of the model, and |
| 816 |
c pl(i,j,np+1) is the surface pressure. |
| 817 |
c (5) the pressure levels ict and icb correspond approximately |
| 818 |
c to 400 and 700 mb. |
| 819 |
c |
| 820 |
c************************************************************************** |
| 821 |
|
| 822 |
implicit none |
| 823 |
|
| 824 |
c-----Explicit Inline Directives |
| 825 |
|
| 826 |
#if CRAY |
| 827 |
#if f77 |
| 828 |
cfpp$ expand (expmn) |
| 829 |
#endif |
| 830 |
#endif |
| 831 |
real expmn |
| 832 |
|
| 833 |
c-----input parameters |
| 834 |
|
| 835 |
integer m,n,ndim,np,ict,icb |
| 836 |
real pl(m,ndim,np+1),ta(m,ndim,np),wa(m,ndim,np),oa(m,ndim,np) |
| 837 |
real taucld(m,ndim,np,2),reff(m,ndim,np,2) |
| 838 |
real fcld(m,ndim,np),taual(m,ndim,np) |
| 839 |
real rsirbm(m,ndim),rsirdf(m,ndim), |
| 840 |
* rsuvbm(m,ndim),rsuvdf(m,ndim),cosz(m,ndim),co2 |
| 841 |
|
| 842 |
c-----output parameters |
| 843 |
|
| 844 |
real flx(m,ndim,np+1),flc(m,ndim,np+1) |
| 845 |
real fdirir(m,ndim),fdifir(m,ndim) |
| 846 |
real fdirpar(m,ndim),fdifpar(m,ndim) |
| 847 |
real fdiruv(m,ndim),fdifuv(m,ndim) |
| 848 |
|
| 849 |
c-----temporary array |
| 850 |
|
| 851 |
integer i,j,k,ik |
| 852 |
real cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np) |
| 853 |
real dp(m,n,np),wh(m,n,np),oh(m,n,np),scal(m,n,np) |
| 854 |
real swh(m,n,np+1),so2(m,n,np+1),df(m,n,np+1) |
| 855 |
real sdf(m,n),sclr(m,n),csm(m,n),taux,x |
| 856 |
|
| 857 |
c----------------------------------------------------------------- |
| 858 |
|
| 859 |
do j= 1, n |
| 860 |
do i= 1, m |
| 861 |
|
| 862 |
swh(i,j,1)=0. |
| 863 |
so2(i,j,1)=0. |
| 864 |
|
| 865 |
c-----csm is the effective secant of the solar zenith angle |
| 866 |
c see equation (12) of Lacis and Hansen (1974, JAS) |
| 867 |
|
| 868 |
csm(i,j)=35./sqrt(1224.*cosz(i,j)*cosz(i,j)+1.) |
| 869 |
|
| 870 |
enddo |
| 871 |
enddo |
| 872 |
|
| 873 |
|
| 874 |
do k= 1, np |
| 875 |
do j= 1, n |
| 876 |
do i= 1, m |
| 877 |
|
| 878 |
c-----compute layer thickness and pressure-scaling function. |
| 879 |
c indices for the surface level and surface layer |
| 880 |
c are np+1 and np, respectively. |
| 881 |
|
| 882 |
dp(i,j,k)=pl(i,j,k+1)-pl(i,j,k) |
| 883 |
scal(i,j,k)=dp(i,j,k)*(.5*(pl(i,j,k)+pl(i,j,k+1))/300.)**.8 |
| 884 |
|
| 885 |
c-----compute scaled water vapor amount, unit is g/cm**2 |
| 886 |
|
| 887 |
wh(i,j,k)=1.02*wa(i,j,k)*scal(i,j,k)* |
| 888 |
* (1.+0.00135*(ta(i,j,k)-240.)) |
| 889 |
swh(i,j,k+1)=swh(i,j,k)+wh(i,j,k) |
| 890 |
|
| 891 |
c-----compute ozone amount, unit is (cm-atm)stp. |
| 892 |
|
| 893 |
oh(i,j,k)=1.02*oa(i,j,k)*dp(i,j,k)*466.7 |
| 894 |
|
| 895 |
enddo |
| 896 |
enddo |
| 897 |
enddo |
| 898 |
|
| 899 |
|
| 900 |
c-----scale cloud optical thickness in each layer from taucld (with |
| 901 |
c cloud amount fcld) to tauclb and tauclf (with cloud amount cc). |
| 902 |
c tauclb is the scaled optical thickness for beam radiation and |
| 903 |
c tauclf is for diffuse radiation. |
| 904 |
|
| 905 |
call cldscale(m,n,ndim,np,cosz,fcld,taucld,ict,icb, |
| 906 |
* cc,tauclb,tauclf) |
| 907 |
|
| 908 |
c-----initialize fluxes for all-sky (flx), clear-sky (flc), and |
| 909 |
c flux reduction (df) |
| 910 |
|
| 911 |
do k=1, np+1 |
| 912 |
do j=1, n |
| 913 |
do i=1, m |
| 914 |
flx(i,j,k)=0. |
| 915 |
flc(i,j,k)=0. |
| 916 |
df(i,j,k)=0. |
| 917 |
enddo |
| 918 |
enddo |
| 919 |
enddo |
| 920 |
|
| 921 |
c-----compute solar ir fluxes |
| 922 |
|
| 923 |
call solir (m,n,ndim,np,wh,taucld,tauclb,tauclf,reff,ict,icb |
| 924 |
* ,fcld,cc,taual,csm,rsirbm,rsirdf,flx,flc,fdirir,fdifir) |
| 925 |
|
| 926 |
c-----compute and update uv and par fluxes |
| 927 |
|
| 928 |
call soluv (m,n,ndim,np,oh,dp,taucld,tauclb,tauclf,reff,ict,icb |
| 929 |
* ,fcld,cc,taual,csm,rsuvbm,rsuvdf,flx,flc |
| 930 |
* ,fdirpar,fdifpar,fdiruv,fdifuv) |
| 931 |
|
| 932 |
c-----compute scaled amount of o2 (so2), unit is (cm-atm)stp. |
| 933 |
|
| 934 |
do k= 1, np |
| 935 |
do j= 1, n |
| 936 |
do i= 1, m |
| 937 |
so2(i,j,k+1)=so2(i,j,k)+165.22*scal(i,j,k) |
| 938 |
enddo |
| 939 |
enddo |
| 940 |
enddo |
| 941 |
|
| 942 |
c-----compute flux reduction due to oxygen following |
| 943 |
c chou (J. climate, 1990). The fraction 0.0287 is the |
| 944 |
c extraterrestrial solar flux in the o2 bands. |
| 945 |
|
| 946 |
do k= 2, np+1 |
| 947 |
do j= 1, n |
| 948 |
do i= 1, m |
| 949 |
x=so2(i,j,k)*csm(i,j) |
| 950 |
df(i,j,k)=df(i,j,k)+0.0287*(1.-expmn(-0.00027*sqrt(x))) |
| 951 |
enddo |
| 952 |
enddo |
| 953 |
enddo |
| 954 |
|
| 955 |
c-----compute scaled amounts for co2 (so2). unit is (cm-atm)stp. |
| 956 |
|
| 957 |
do k= 1, np |
| 958 |
do j= 1, n |
| 959 |
do i= 1, m |
| 960 |
so2(i,j,k+1)=so2(i,j,k)+co2*789.*scal(i,j,k) |
| 961 |
enddo |
| 962 |
enddo |
| 963 |
enddo |
| 964 |
|
| 965 |
c-----compute and update flux reduction due to co2 following |
| 966 |
c chou (J. Climate, 1990) |
| 967 |
|
| 968 |
call flxco2(m,n,np,so2,swh,csm,df) |
| 969 |
|
| 970 |
c-----adjust for the effect of o2 cnd co2 on clear-sky fluxes. |
| 971 |
|
| 972 |
do k= 2, np+1 |
| 973 |
do j= 1, n |
| 974 |
do i= 1, m |
| 975 |
flc(i,j,k)=flc(i,j,k)-df(i,j,k) |
| 976 |
enddo |
| 977 |
enddo |
| 978 |
enddo |
| 979 |
|
| 980 |
c-----adjust for the all-sky fluxes due to o2 and co2. It is |
| 981 |
c assumed that o2 and co2 have no effects on solar radiation |
| 982 |
c below clouds. clouds are assumed randomly overlapped. |
| 983 |
|
| 984 |
do j=1,n |
| 985 |
do i=1,m |
| 986 |
sdf(i,j)=0.0 |
| 987 |
sclr(i,j)=1.0 |
| 988 |
enddo |
| 989 |
enddo |
| 990 |
|
| 991 |
do k=1,np |
| 992 |
do j=1,n |
| 993 |
do i=1,m |
| 994 |
|
| 995 |
c-----sclr is the fraction of clear sky. |
| 996 |
c sdf is the flux reduction below clouds. |
| 997 |
|
| 998 |
if(fcld(i,j,k).gt.0.01) then |
| 999 |
sdf(i,j)=sdf(i,j)+df(i,j,k)*sclr(i,j)*fcld(i,j,k) |
| 1000 |
sclr(i,j)=sclr(i,j)*(1.-fcld(i,j,k)) |
| 1001 |
endif |
| 1002 |
flx(i,j,k+1)=flx(i,j,k+1)-sdf(i,j)-df(i,j,k+1)*sclr(i,j) |
| 1003 |
|
| 1004 |
enddo |
| 1005 |
enddo |
| 1006 |
enddo |
| 1007 |
|
| 1008 |
c-----adjust for the direct downward ir flux. |
| 1009 |
do j= 1, n |
| 1010 |
do i= 1, m |
| 1011 |
x = (1.-rsirdf(i,j))*fdifir(i,j) + (1.-rsirbm(i,j))*fdirir(i,j) |
| 1012 |
x = (-sdf(i,j)-df(i,j,np+1)*sclr(i,j))/x |
| 1013 |
fdirir(i,j)=fdirir(i,j)*(1.+x) |
| 1014 |
fdifir(i,j)=fdifir(i,j)*(1.+x) |
| 1015 |
enddo |
| 1016 |
enddo |
| 1017 |
|
| 1018 |
return |
| 1019 |
end |
| 1020 |
|
| 1021 |
c******************************************************************** |
| 1022 |
|
| 1023 |
subroutine cldscale (m,n,ndim,np,cosz,fcld,taucld,ict,icb, |
| 1024 |
* cc,tauclb,tauclf) |
| 1025 |
|
| 1026 |
c******************************************************************** |
| 1027 |
c |
| 1028 |
c This subroutine computes the covers of high, middle, and |
| 1029 |
c low clouds and scales the cloud optical thickness. |
| 1030 |
c |
| 1031 |
c To simplify calculations in a cloudy atmosphere, clouds are |
| 1032 |
c grouped into high, middle and low clouds separated by the levels |
| 1033 |
c ict and icb (level 1 is the top of the atmosphere). |
| 1034 |
c |
| 1035 |
c Within each of the three groups, clouds are assumed maximally |
| 1036 |
c overlapped, and the cloud cover (cc) of a group is the maximum |
| 1037 |
c cloud cover of all the layers in the group. The optical thickness |
| 1038 |
c (taucld) of a given layer is then scaled to new values (tauclb and |
| 1039 |
c tauclf) so that the layer reflectance corresponding to the cloud |
| 1040 |
c cover cc is the same as the original reflectance with optical |
| 1041 |
c thickness taucld and cloud cover fcld. |
| 1042 |
c |
| 1043 |
c---input parameters |
| 1044 |
c |
| 1045 |
c number of grid intervals in zonal direction (m) |
| 1046 |
c number of grid intervals in meridional direction (n) |
| 1047 |
c maximum number of grid intervals in meridional direction (ndim) |
| 1048 |
c number of atmospheric layers (np) |
| 1049 |
c cosine of the solar zenith angle (cosz) |
| 1050 |
c fractional cloud cover (fcld) |
| 1051 |
c cloud optical thickness (taucld) |
| 1052 |
c index separating high and middle clouds (ict) |
| 1053 |
c index separating middle and low clouds (icb) |
| 1054 |
c |
| 1055 |
c---output parameters |
| 1056 |
c |
| 1057 |
c fractional cover of high, middle, and low clouds (cc) |
| 1058 |
c scaled cloud optical thickness for beam radiation (tauclb) |
| 1059 |
c scaled cloud optical thickness for diffuse radiation (tauclf) |
| 1060 |
c |
| 1061 |
c******************************************************************** |
| 1062 |
|
| 1063 |
implicit none |
| 1064 |
|
| 1065 |
c-----input parameters |
| 1066 |
|
| 1067 |
integer m,n,ndim,np,ict,icb |
| 1068 |
real cosz(m,ndim),fcld(m,ndim,np),taucld(m,ndim,np,2) |
| 1069 |
|
| 1070 |
c-----output parameters |
| 1071 |
|
| 1072 |
real cc(m,n,3),tauclb(m,n,np),tauclf(m,n,np) |
| 1073 |
|
| 1074 |
c-----temporary variables |
| 1075 |
|
| 1076 |
integer i,j,k,im,it,ia,kk |
| 1077 |
real fm,ft,fa,xai,taucl,taux |
| 1078 |
|
| 1079 |
c-----pre-computed table |
| 1080 |
|
| 1081 |
integer nm,nt,na |
| 1082 |
parameter (nm=11,nt=9,na=11) |
| 1083 |
real dm,dt,da,t1,caib(nm,nt,na),caif(nt,na) |
| 1084 |
parameter (dm=0.1,dt=0.30103,da=0.1,t1=-0.9031) |
| 1085 |
|
| 1086 |
c-----include the pre-computed table for cai |
| 1087 |
|
| 1088 |
include 'cai.dat' |
| 1089 |
save caib,caif |
| 1090 |
|
| 1091 |
|
| 1092 |
c-----clouds within each of the high, middle, and low clouds are |
| 1093 |
c assumed maximally overlapped, and the cloud cover (cc) |
| 1094 |
c for a group is the maximum cloud cover of all the layers |
| 1095 |
c in the group |
| 1096 |
|
| 1097 |
do j=1,n |
| 1098 |
do i=1,m |
| 1099 |
cc(i,j,1)=0.0 |
| 1100 |
cc(i,j,2)=0.0 |
| 1101 |
cc(i,j,3)=0.0 |
| 1102 |
enddo |
| 1103 |
enddo |
| 1104 |
|
| 1105 |
do k=1,ict-1 |
| 1106 |
do j=1,n |
| 1107 |
do i=1,m |
| 1108 |
cc(i,j,1)=max(cc(i,j,1),fcld(i,j,k)) |
| 1109 |
enddo |
| 1110 |
enddo |
| 1111 |
enddo |
| 1112 |
|
| 1113 |
do k=ict,icb-1 |
| 1114 |
do j=1,n |
| 1115 |
do i=1,m |
| 1116 |
cc(i,j,2)=max(cc(i,j,2),fcld(i,j,k)) |
| 1117 |
enddo |
| 1118 |
enddo |
| 1119 |
enddo |
| 1120 |
|
| 1121 |
do k=icb,np |
| 1122 |
do j=1,n |
| 1123 |
do i=1,m |
| 1124 |
cc(i,j,3)=max(cc(i,j,3),fcld(i,j,k)) |
| 1125 |
enddo |
| 1126 |
enddo |
| 1127 |
enddo |
| 1128 |
|
| 1129 |
c-----scale the cloud optical thickness. |
| 1130 |
c taucld(i,j,k,1) is the optical thickness for ice particles, and |
| 1131 |
c taucld(i,j,k,2) is the optical thickness for liquid particles. |
| 1132 |
|
| 1133 |
do k=1,np |
| 1134 |
|
| 1135 |
if(k.lt.ict) then |
| 1136 |
kk=1 |
| 1137 |
elseif(k.ge.ict .and. k.lt.icb) then |
| 1138 |
kk=2 |
| 1139 |
else |
| 1140 |
kk=3 |
| 1141 |
endif |
| 1142 |
|
| 1143 |
do j=1,n |
| 1144 |
do i=1,m |
| 1145 |
|
| 1146 |
tauclb(i,j,k) = 0.0 |
| 1147 |
tauclf(i,j,k) = 0.0 |
| 1148 |
|
| 1149 |
taux=taucld(i,j,k,1)+taucld(i,j,k,2) |
| 1150 |
if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then |
| 1151 |
|
| 1152 |
c-----normalize cloud cover |
| 1153 |
|
| 1154 |
fa=fcld(i,j,k)/cc(i,j,kk) |
| 1155 |
|
| 1156 |
c-----table look-up |
| 1157 |
|
| 1158 |
taux=min(taux,32.) |
| 1159 |
|
| 1160 |
fm=cosz(i,j)/dm |
| 1161 |
ft=(log10(taux)-t1)/dt |
| 1162 |
fa=fa/da |
| 1163 |
|
| 1164 |
im=int(fm+1.5) |
| 1165 |
it=int(ft+1.5) |
| 1166 |
ia=int(fa+1.5) |
| 1167 |
|
| 1168 |
im=max(im,2) |
| 1169 |
it=max(it,2) |
| 1170 |
ia=max(ia,2) |
| 1171 |
|
| 1172 |
im=min(im,nm-1) |
| 1173 |
it=min(it,nt-1) |
| 1174 |
ia=min(ia,na-1) |
| 1175 |
|
| 1176 |
fm=fm-float(im-1) |
| 1177 |
ft=ft-float(it-1) |
| 1178 |
fa=fa-float(ia-1) |
| 1179 |
|
| 1180 |
c-----scale cloud optical thickness for beam radiation. |
| 1181 |
c the scaling factor, xai, is a function of the solar zenith |
| 1182 |
c angle, optical thickness, and cloud cover. |
| 1183 |
|
| 1184 |
xai= (-caib(im-1,it,ia)*(1.-fm)+ |
| 1185 |
* caib(im+1,it,ia)*(1.+fm))*fm*.5+caib(im,it,ia)*(1.-fm*fm) |
| 1186 |
|
| 1187 |
xai=xai+(-caib(im,it-1,ia)*(1.-ft)+ |
| 1188 |
* caib(im,it+1,ia)*(1.+ft))*ft*.5+caib(im,it,ia)*(1.-ft*ft) |
| 1189 |
|
| 1190 |
xai=xai+(-caib(im,it,ia-1)*(1.-fa)+ |
| 1191 |
* caib(im,it,ia+1)*(1.+fa))*fa*.5+caib(im,it,ia)*(1.-fa*fa) |
| 1192 |
|
| 1193 |
xai= xai-2.*caib(im,it,ia) |
| 1194 |
xai=max(xai,0.0) |
| 1195 |
|
| 1196 |
tauclb(i,j,k) = taux*xai |
| 1197 |
|
| 1198 |
c-----scale cloud optical thickness for diffuse radiation. |
| 1199 |
c the scaling factor, xai, is a function of the cloud optical |
| 1200 |
c thickness and cover but not the solar zenith angle. |
| 1201 |
|
| 1202 |
xai= (-caif(it-1,ia)*(1.-ft)+ |
| 1203 |
* caif(it+1,ia)*(1.+ft))*ft*.5+caif(it,ia)*(1.-ft*ft) |
| 1204 |
|
| 1205 |
xai=xai+(-caif(it,ia-1)*(1.-fa)+ |
| 1206 |
* caif(it,ia+1)*(1.+fa))*fa*.5+caif(it,ia)*(1.-fa*fa) |
| 1207 |
|
| 1208 |
xai= xai-caif(it,ia) |
| 1209 |
xai=max(xai,0.0) |
| 1210 |
|
| 1211 |
tauclf(i,j,k) = taux*xai |
| 1212 |
|
| 1213 |
endif |
| 1214 |
|
| 1215 |
enddo |
| 1216 |
enddo |
| 1217 |
enddo |
| 1218 |
|
| 1219 |
return |
| 1220 |
end |
| 1221 |
c*********************************************************************** |
| 1222 |
|
| 1223 |
subroutine solir (m,n,ndim,np,wh,taucld,tauclb,tauclf,reff, |
| 1224 |
* ict,icb,fcld,cc,taual,csm,rsirbm,rsirdf,flx,flc,fdirir,fdifir) |
| 1225 |
|
| 1226 |
c************************************************************************ |
| 1227 |
c compute solar flux in the infrared region. The spectrum is divided |
| 1228 |
c into three bands: |
| 1229 |
c |
| 1230 |
c band wavenumber(/cm) wavelength (micron) |
| 1231 |
c 1 1000-4400 2.27-10.0 |
| 1232 |
c 2 4400-8200 1.22-2.27 |
| 1233 |
c 3 8200-14300 0.70-1.22 |
| 1234 |
c |
| 1235 |
c----- Input parameters: units size |
| 1236 |
c |
| 1237 |
c number of soundings in zonal direction (m) n/d 1 |
| 1238 |
c number of soundings in meridional direction (n) n/d 1 |
| 1239 |
c maximum number of soundings in n/d 1 |
| 1240 |
c meridional direction (ndim) |
| 1241 |
c number of atmospheric layers (np) n/d 1 |
| 1242 |
c layer water vapor content (wh) gm/cm^2 m*n*np |
| 1243 |
c cloud optical thickness (taucld) n/d m*ndim*np*2 |
| 1244 |
c index 1 for ice paticles |
| 1245 |
c index 2 for liquid particles |
| 1246 |
c scaled cloud optical thickness n/d m*n*np |
| 1247 |
c for beam radiation (tauclb) |
| 1248 |
c scaled cloud optical thickness n/d m*n*np |
| 1249 |
c for diffuse radiation (tauclf) |
| 1250 |
c effective cloud-particle size (reff) micrometer m*ndim*np*2 |
| 1251 |
c index 1 for ice paticles |
| 1252 |
c index 2 for liquid particles |
| 1253 |
c level index separating high and n/d m*n |
| 1254 |
c middle clouds (ict) |
| 1255 |
c level index separating middle and n/d m*n |
| 1256 |
c low clouds (icb) |
| 1257 |
c cloud amount (fcld) fraction m*ndim*np |
| 1258 |
c cloud amount of high, middle, and n/d m*n*3 |
| 1259 |
c low clouds (cc) |
| 1260 |
c aerosol optical thickness (taual) n/d m*ndim*np |
| 1261 |
c cosecant of the solar zenith angle (csm) n/d m*n |
| 1262 |
c near ir surface albedo for beam fraction m*ndim |
| 1263 |
c radiation (rsirbm) |
| 1264 |
c near ir surface albedo for diffuse fraction m*ndim |
| 1265 |
c radiation (rsirdf) |
| 1266 |
c |
| 1267 |
c----- output (updated) parameters: |
| 1268 |
c |
| 1269 |
c all-sky flux (downward-upward) (flx) fraction m*ndim*(np+1) |
| 1270 |
c clear-sky flux (downward-upward) (flc) fraction m*ndim*(np+1) |
| 1271 |
c all-sky direct downward ir flux at |
| 1272 |
c the surface (fdirir) fraction m*ndim |
| 1273 |
c all-sky diffuse downward ir flux at |
| 1274 |
c the surface (fdifir) fraction m*ndim |
| 1275 |
c |
| 1276 |
c----- note: the following parameters must be specified by users: |
| 1277 |
c aerosol single scattering albedo (ssaal) n/d nband |
| 1278 |
c aerosol asymmetry factor (asyal) n/d nband |
| 1279 |
c |
| 1280 |
c************************************************************************* |
| 1281 |
|
| 1282 |
implicit none |
| 1283 |
|
| 1284 |
c-----Explicit Inline Directives |
| 1285 |
|
| 1286 |
#if CRAY |
| 1287 |
#if f77 |
| 1288 |
cfpp$ expand (deledd) |
| 1289 |
cfpp$ expand (sagpol) |
| 1290 |
cfpp$ expand (expmn) |
| 1291 |
#endif |
| 1292 |
#endif |
| 1293 |
real expmn |
| 1294 |
|
| 1295 |
c-----input parameters |
| 1296 |
|
| 1297 |
integer m,n,ndim,np,ict,icb |
| 1298 |
real taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np) |
| 1299 |
real tauclb(m,n,np),tauclf(m,n,np),cc(m,n,3) |
| 1300 |
real rsirbm(m,ndim),rsirdf(m,ndim) |
| 1301 |
real wh(m,n,np),taual(m,ndim,np),csm(m,n) |
| 1302 |
|
| 1303 |
c-----output (updated) parameters |
| 1304 |
|
| 1305 |
real flx(m,ndim,np+1),flc(m,ndim,np+1) |
| 1306 |
real fdirir(m,ndim),fdifir(m,ndim) |
| 1307 |
|
| 1308 |
c-----static parameters |
| 1309 |
|
| 1310 |
integer nk,nband |
| 1311 |
parameter (nk=10,nband=3) |
| 1312 |
real xk(nk),hk(nband,nk),ssaal(nband),asyal(nband) |
| 1313 |
real aia(nband,3),awa(nband,3),aig(nband,3),awg(nband,3) |
| 1314 |
|
| 1315 |
c-----temporary array |
| 1316 |
|
| 1317 |
integer ib,ik,i,j,k |
| 1318 |
real ssacl(m,n,np),asycl(m,n,np) |
| 1319 |
real rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2), |
| 1320 |
* rs(m,n,np+1,2),ts(m,n,np+1,2) |
| 1321 |
real rssab(m,n,np+1),rabx(m,n,np+1),rsabx(m,n,np+1) |
| 1322 |
real fall(m,n,np+1),fclr(m,n,np+1) |
| 1323 |
real fsdir(m,n),fsdif(m,n) |
| 1324 |
|
| 1325 |
real tauwv,tausto,ssatau,asysto,tauto,ssato,asyto |
| 1326 |
real taux,reff1,reff2,w1,w2,g1,g2 |
| 1327 |
real ssaclt(m,n),asyclt(m,n) |
| 1328 |
real rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n) |
| 1329 |
real rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n) |
| 1330 |
|
| 1331 |
c-----water vapor absorption coefficient for 10 k-intervals. |
| 1332 |
c unit: cm^2/gm |
| 1333 |
|
| 1334 |
data xk/ |
| 1335 |
1 0.0010, 0.0133, 0.0422, 0.1334, 0.4217, |
| 1336 |
2 1.334, 5.623, 31.62, 177.8, 1000.0/ |
| 1337 |
|
| 1338 |
c-----water vapor k-distribution function, |
| 1339 |
c the sum of hk is 0.52926. unit: fraction |
| 1340 |
|
| 1341 |
data hk/ |
| 1342 |
1 .01074,.08236,.20673, .00360,.01157,.03497, |
| 1343 |
2 .00411,.01133,.03011, .00421,.01143,.02260, |
| 1344 |
3 .00389,.01240,.01336, .00326,.01258,.00696, |
| 1345 |
4 .00499,.01381,.00441, .00465,.00650,.00115, |
| 1346 |
5 .00245,.00244,.00026, .00145,.00094,.00000/ |
| 1347 |
|
| 1348 |
c-----aerosol single-scattering albedo and asymmetry factor |
| 1349 |
|
| 1350 |
data ssaal,asyal/nband*0.999,nband*0.850/ |
| 1351 |
|
| 1352 |
c-----coefficients for computing the single scattering albedo of |
| 1353 |
c ice clouds from ssa=1-(aia(*,1)+aia(*,2)*reff+aia(*,3)*reff**2) |
| 1354 |
|
| 1355 |
data aia/ |
| 1356 |
1 .08938331, .00215346,-.00000260, |
| 1357 |
2 .00299387, .00073709, .00000746, |
| 1358 |
3 -.00001038,-.00000134, .00000000/ |
| 1359 |
|
| 1360 |
c-----coefficients for computing the single scattering albedo of |
| 1361 |
c liquid clouds from ssa=1-(awa(*,1)+awa(*,2)*reff+awa(*,3)*reff**2) |
| 1362 |
|
| 1363 |
data awa/ |
| 1364 |
1 .01209318,-.00019934, .00000007, |
| 1365 |
2 .01784739, .00088757, .00000845, |
| 1366 |
3 -.00036910,-.00000650,-.00000004/ |
| 1367 |
|
| 1368 |
c-----coefficients for computing the asymmetry factor of ice clouds |
| 1369 |
c from asycl=aig(*,1)+aig(*,2)*reff+aig(*,3)*reff**2 |
| 1370 |
|
| 1371 |
data aig/ |
| 1372 |
1 .84090400, .76098937, .74935228, |
| 1373 |
2 .00126222, .00141864, .00119715, |
| 1374 |
3 -.00000385,-.00000396,-.00000367/ |
| 1375 |
|
| 1376 |
c-----coefficients for computing the asymmetry factor of liquid clouds |
| 1377 |
c from asycl=awg(*,1)+awg(*,2)*reff+awg(*,3)*reff**2 |
| 1378 |
|
| 1379 |
data awg/ |
| 1380 |
1 .83530748, .74513197, .79375035, |
| 1381 |
2 .00257181, .01370071, .00832441, |
| 1382 |
3 .00005519,-.00038203,-.00023263/ |
| 1383 |
|
| 1384 |
c-----initialize surface fluxes, reflectances, and transmittances |
| 1385 |
|
| 1386 |
do j= 1, n |
| 1387 |
do i= 1, m |
| 1388 |
fdirir(i,j)=0.0 |
| 1389 |
fdifir(i,j)=0.0 |
| 1390 |
rr(i,j,np+1,1)=rsirbm(i,j) |
| 1391 |
rr(i,j,np+1,2)=rsirbm(i,j) |
| 1392 |
rs(i,j,np+1,1)=rsirdf(i,j) |
| 1393 |
rs(i,j,np+1,2)=rsirdf(i,j) |
| 1394 |
td(i,j,np+1,1)=0.0 |
| 1395 |
td(i,j,np+1,2)=0.0 |
| 1396 |
tt(i,j,np+1,1)=0.0 |
| 1397 |
tt(i,j,np+1,2)=0.0 |
| 1398 |
ts(i,j,np+1,1)=0.0 |
| 1399 |
ts(i,j,np+1,2)=0.0 |
| 1400 |
enddo |
| 1401 |
enddo |
| 1402 |
|
| 1403 |
c-----integration over spectral bands |
| 1404 |
|
| 1405 |
do 100 ib=1,nband |
| 1406 |
|
| 1407 |
c-----compute cloud single scattering albedo and asymmetry factor |
| 1408 |
c for a mixture of ice and liquid particles. |
| 1409 |
|
| 1410 |
do k= 1, np |
| 1411 |
|
| 1412 |
do j= 1, n |
| 1413 |
do i= 1, m |
| 1414 |
|
| 1415 |
ssaclt(i,j)=1.0 |
| 1416 |
asyclt(i,j)=1.0 |
| 1417 |
|
| 1418 |
taux=taucld(i,j,k,1)+taucld(i,j,k,2) |
| 1419 |
if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then |
| 1420 |
|
| 1421 |
reff1=min(reff(i,j,k,1),130.) |
| 1422 |
reff2=min(reff(i,j,k,2),20.0) |
| 1423 |
|
| 1424 |
w1=(1.-(aia(ib,1)+(aia(ib,2)+ |
| 1425 |
* aia(ib,3)*reff1)*reff1))*taucld(i,j,k,1) |
| 1426 |
w2=(1.-(awa(ib,1)+(awa(ib,2)+ |
| 1427 |
* awa(ib,3)*reff2)*reff2))*taucld(i,j,k,2) |
| 1428 |
ssaclt(i,j)=(w1+w2)/taux |
| 1429 |
|
| 1430 |
g1=(aig(ib,1)+(aig(ib,2)+aig(ib,3)*reff1)*reff1)*w1 |
| 1431 |
g2=(awg(ib,1)+(awg(ib,2)+awg(ib,3)*reff2)*reff2)*w2 |
| 1432 |
asyclt(i,j)=(g1+g2)/(w1+w2) |
| 1433 |
|
| 1434 |
endif |
| 1435 |
|
| 1436 |
enddo |
| 1437 |
enddo |
| 1438 |
|
| 1439 |
do j=1,n |
| 1440 |
do i=1,m |
| 1441 |
ssacl(i,j,k)=ssaclt(i,j) |
| 1442 |
enddo |
| 1443 |
enddo |
| 1444 |
do j=1,n |
| 1445 |
do i=1,m |
| 1446 |
asycl(i,j,k)=asyclt(i,j) |
| 1447 |
enddo |
| 1448 |
enddo |
| 1449 |
|
| 1450 |
enddo |
| 1451 |
|
| 1452 |
c-----integration over the k-distribution function |
| 1453 |
|
| 1454 |
do 200 ik=1,nk |
| 1455 |
|
| 1456 |
do 300 k= 1, np |
| 1457 |
|
| 1458 |
do j= 1, n |
| 1459 |
do i= 1, m |
| 1460 |
|
| 1461 |
tauwv=xk(ik)*wh(i,j,k) |
| 1462 |
|
| 1463 |
c-----compute total optical thickness, single scattering albedo, |
| 1464 |
c and asymmetry factor. |
| 1465 |
|
| 1466 |
tausto=tauwv+taual(i,j,k)+1.0e-8 |
| 1467 |
ssatau=ssaal(ib)*taual(i,j,k) |
| 1468 |
asysto=asyal(ib)*ssaal(ib)*taual(i,j,k) |
| 1469 |
|
| 1470 |
c-----compute reflectance and transmittance for cloudless layers |
| 1471 |
|
| 1472 |
tauto=tausto |
| 1473 |
ssato=ssatau/tauto+1.0e-8 |
| 1474 |
|
| 1475 |
c if (ssato .gt. 0.001) then |
| 1476 |
|
| 1477 |
c ssato=min(ssato,0.999999) |
| 1478 |
c asyto=asysto/(ssato*tauto) |
| 1479 |
|
| 1480 |
c call deledd(tauto,ssato,asyto,csm(i,j), |
| 1481 |
c * rr1t(i,j),tt1t(i,j),td1t(i,j)) |
| 1482 |
|
| 1483 |
c call sagpol (tauto,ssato,asyto,rs1t(i,j),ts1t(i,j)) |
| 1484 |
|
| 1485 |
c else |
| 1486 |
|
| 1487 |
td1t(i,j)=expmn(-tauto*csm(i,j)) |
| 1488 |
ts1t(i,j)=expmn(-1.66*tauto) |
| 1489 |
tt1t(i,j)=0.0 |
| 1490 |
rr1t(i,j)=0.0 |
| 1491 |
rs1t(i,j)=0.0 |
| 1492 |
|
| 1493 |
c endif |
| 1494 |
|
| 1495 |
c-----compute reflectance and transmittance for cloud layers |
| 1496 |
|
| 1497 |
if (tauclb(i,j,k) .lt. 0.01) then |
| 1498 |
|
| 1499 |
rr2t(i,j)=rr1t(i,j) |
| 1500 |
tt2t(i,j)=tt1t(i,j) |
| 1501 |
td2t(i,j)=td1t(i,j) |
| 1502 |
rs2t(i,j)=rs1t(i,j) |
| 1503 |
ts2t(i,j)=ts1t(i,j) |
| 1504 |
|
| 1505 |
else |
| 1506 |
|
| 1507 |
tauto=tausto+tauclb(i,j,k) |
| 1508 |
ssato=(ssatau+ssacl(i,j,k)*tauclb(i,j,k))/tauto+1.0e-8 |
| 1509 |
ssato=min(ssato,0.999999) |
| 1510 |
asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclb(i,j,k))/ |
| 1511 |
* (ssato*tauto) |
| 1512 |
|
| 1513 |
call deledd(tauto,ssato,asyto,csm(i,j), |
| 1514 |
* rr2t(i,j),tt2t(i,j),td2t(i,j)) |
| 1515 |
|
| 1516 |
tauto=tausto+tauclf(i,j,k) |
| 1517 |
ssato=(ssatau+ssacl(i,j,k)*tauclf(i,j,k))/tauto+1.0e-8 |
| 1518 |
ssato=min(ssato,0.999999) |
| 1519 |
asyto=(asysto+asycl(i,j,k)*ssacl(i,j,k)*tauclf(i,j,k))/ |
| 1520 |
* (ssato*tauto) |
| 1521 |
|
| 1522 |
call sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j)) |
| 1523 |
|
| 1524 |
endif |
| 1525 |
|
| 1526 |
enddo |
| 1527 |
enddo |
| 1528 |
|
| 1529 |
c-----the following code appears to be awkward, but it is efficient |
| 1530 |
c in certain parallel processors. |
| 1531 |
|
| 1532 |
do j=1,n |
| 1533 |
do i=1,m |
| 1534 |
rr(i,j,k,1)=rr1t(i,j) |
| 1535 |
enddo |
| 1536 |
enddo |
| 1537 |
do j=1,n |
| 1538 |
do i=1,m |
| 1539 |
tt(i,j,k,1)=tt1t(i,j) |
| 1540 |
enddo |
| 1541 |
enddo |
| 1542 |
do j=1,n |
| 1543 |
do i=1,m |
| 1544 |
td(i,j,k,1)=td1t(i,j) |
| 1545 |
enddo |
| 1546 |
enddo |
| 1547 |
do j=1,n |
| 1548 |
do i=1,m |
| 1549 |
rs(i,j,k,1)=rs1t(i,j) |
| 1550 |
enddo |
| 1551 |
enddo |
| 1552 |
do j=1,n |
| 1553 |
do i=1,m |
| 1554 |
ts(i,j,k,1)=ts1t(i,j) |
| 1555 |
enddo |
| 1556 |
enddo |
| 1557 |
|
| 1558 |
do j=1,n |
| 1559 |
do i=1,m |
| 1560 |
rr(i,j,k,2)=rr2t(i,j) |
| 1561 |
enddo |
| 1562 |
enddo |
| 1563 |
do j=1,n |
| 1564 |
do i=1,m |
| 1565 |
tt(i,j,k,2)=tt2t(i,j) |
| 1566 |
enddo |
| 1567 |
enddo |
| 1568 |
do j=1,n |
| 1569 |
do i=1,m |
| 1570 |
td(i,j,k,2)=td2t(i,j) |
| 1571 |
enddo |
| 1572 |
enddo |
| 1573 |
do j=1,n |
| 1574 |
do i=1,m |
| 1575 |
rs(i,j,k,2)=rs2t(i,j) |
| 1576 |
enddo |
| 1577 |
enddo |
| 1578 |
do j=1,n |
| 1579 |
do i=1,m |
| 1580 |
ts(i,j,k,2)=ts2t(i,j) |
| 1581 |
enddo |
| 1582 |
enddo |
| 1583 |
|
| 1584 |
300 continue |
| 1585 |
|
| 1586 |
c-----flux calculations |
| 1587 |
|
| 1588 |
call cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts, |
| 1589 |
* fclr,fall,fsdir,fsdif) |
| 1590 |
|
| 1591 |
do k= 1, np+1 |
| 1592 |
do j= 1, n |
| 1593 |
do i= 1, m |
| 1594 |
flx(i,j,k) = flx(i,j,k)+fall(i,j,k)*hk(ib,ik) |
| 1595 |
enddo |
| 1596 |
enddo |
| 1597 |
do j= 1, n |
| 1598 |
do i= 1, m |
| 1599 |
flc(i,j,k) = flc(i,j,k)+fclr(i,j,k)*hk(ib,ik) |
| 1600 |
enddo |
| 1601 |
enddo |
| 1602 |
enddo |
| 1603 |
|
| 1604 |
do j= 1, n |
| 1605 |
do i= 1, m |
| 1606 |
fdirir(i,j) = fdirir(i,j)+fsdir(i,j)*hk(ib,ik) |
| 1607 |
fdifir(i,j) = fdifir(i,j)+fsdif(i,j)*hk(ib,ik) |
| 1608 |
enddo |
| 1609 |
enddo |
| 1610 |
|
| 1611 |
200 continue |
| 1612 |
100 continue |
| 1613 |
|
| 1614 |
return |
| 1615 |
end |
| 1616 |
|
| 1617 |
c************************************************************************ |
| 1618 |
|
| 1619 |
subroutine soluv (m,n,ndim,np,oh,dp,taucld,tauclb,tauclf,reff, |
| 1620 |
* ict,icb,fcld,cc,taual,csm,rsuvbm,rsuvdf,flx,flc |
| 1621 |
* ,fdirpar,fdifpar,fdiruv,fdifuv) |
| 1622 |
|
| 1623 |
c************************************************************************ |
| 1624 |
c compute solar fluxes in the uv+visible region. the spectrum is |
| 1625 |
c grouped into 8 bands: |
| 1626 |
c |
| 1627 |
c Band Micrometer |
| 1628 |
c |
| 1629 |
c UV-C 1. .175 - .225 |
| 1630 |
c 2. .225 - .245 |
| 1631 |
c .260 - .280 |
| 1632 |
c 3. .245 - .260 |
| 1633 |
c |
| 1634 |
c UV-B 4. .280 - .295 |
| 1635 |
c 5. .295 - .310 |
| 1636 |
c 6. .310 - .320 |
| 1637 |
c |
| 1638 |
c UV-A 7. .320 - .400 |
| 1639 |
c |
| 1640 |
c PAR 8. .400 - .700 |
| 1641 |
c |
| 1642 |
c----- Input parameters: units size |
| 1643 |
c |
| 1644 |
c number of soundings in zonal direction (m) n/d 1 |
| 1645 |
c number of soundings in meridional direction (n) n/d 1 |
| 1646 |
c maximum number of soundings in n/d 1 |
| 1647 |
c meridional direction (ndim) |
| 1648 |
c number of atmospheric layers (np) n/d 1 |
| 1649 |
c layer ozone content (oh) (cm-atm)stp m*n*np |
| 1650 |
c layer pressure thickness (dp) mb m*n*np |
| 1651 |
c cloud optical thickness (taucld) n/d m*ndim*np*2 |
| 1652 |
c index 1 for ice paticles |
| 1653 |
c index 2 for liquid particles |
| 1654 |
c scaled cloud optical thickness n/d m*n*np |
| 1655 |
c for beam radiation (tauclb) |
| 1656 |
c scaled cloud optical thickness n/d m*n*np |
| 1657 |
c for diffuse radiation (tauclf) |
| 1658 |
c effective cloud-particle size (reff) micrometer m*ndim*np*2 |
| 1659 |
c index 1 for ice paticles |
| 1660 |
c index 2 for liquid particles |
| 1661 |
c level indiex separating high and n/d m*n |
| 1662 |
c middle clouds (ict) |
| 1663 |
c level indiex separating middle and n/d m*n |
| 1664 |
c low clouds (icb) |
| 1665 |
c cloud amount (fcld) fraction m*ndim*np |
| 1666 |
c cloud amounts of high, middle, and n/d m*n*3 |
| 1667 |
c low clouds (cc) |
| 1668 |
c aerosol optical thickness (taual) n/d m*ndim*np |
| 1669 |
c cosecant of the solar zenith angle (csm) n/d m*n |
| 1670 |
c uv+par surface albedo for beam fraction m*ndim |
| 1671 |
c radiation (rsuvbm) |
| 1672 |
c uv+par surface albedo for diffuse fraction m*ndim |
| 1673 |
c radiation (rsuvdf) |
| 1674 |
c |
| 1675 |
c----- output (updated) parameters: |
| 1676 |
c |
| 1677 |
c all-sky net downward flux (flx) fraction m*ndim*(np+1) |
| 1678 |
c clear-sky net downward flux (flc) fraction m*ndim*(np+1) |
| 1679 |
c all-sky direct downward par flux at |
| 1680 |
c the surface (fdirpar) fraction m*ndim |
| 1681 |
c all-sky diffuse downward par flux at |
| 1682 |
c the surface (fdifpar) fraction m*ndim |
| 1683 |
c |
| 1684 |
c----- note: the following parameters must be specified by users: |
| 1685 |
c |
| 1686 |
c aerosol single scattering albedo (ssaal) n/d 1 |
| 1687 |
c aerosol asymmetry factor (asyal) n/d 1 |
| 1688 |
c |
| 1689 |
* |
| 1690 |
c*********************************************************************** |
| 1691 |
|
| 1692 |
implicit none |
| 1693 |
|
| 1694 |
c-----Explicit Inline Directives |
| 1695 |
|
| 1696 |
#if CRAY |
| 1697 |
#if f77 |
| 1698 |
cfpp$ expand (deledd) |
| 1699 |
cfpp$ expand (sagpol) |
| 1700 |
#endif |
| 1701 |
#endif |
| 1702 |
|
| 1703 |
c-----input parameters |
| 1704 |
|
| 1705 |
integer m,n,ndim,np,ict,icb |
| 1706 |
real taucld(m,ndim,np,2),reff(m,ndim,np,2),fcld(m,ndim,np) |
| 1707 |
real tauclb(m,n,np),tauclf(m,n,np),cc(m,n,3) |
| 1708 |
real oh(m,n,np),dp(m,n,np),taual(m,ndim,np) |
| 1709 |
real rsuvbm(m,ndim),rsuvdf(m,ndim),csm(m,n) |
| 1710 |
|
| 1711 |
c-----output (updated) parameter |
| 1712 |
|
| 1713 |
real flx(m,ndim,np+1),flc(m,ndim,np+1) |
| 1714 |
real fdirpar(m,ndim),fdifpar(m,ndim) |
| 1715 |
real fdiruv(m,ndim),fdifuv(m,ndim) |
| 1716 |
|
| 1717 |
c-----static parameters |
| 1718 |
|
| 1719 |
integer nband |
| 1720 |
parameter (nband=8) |
| 1721 |
real hk(nband),xk(nband),ry(nband) |
| 1722 |
real asyal(nband),ssaal(nband),aig(3),awg(3) |
| 1723 |
|
| 1724 |
c-----temporary array |
| 1725 |
|
| 1726 |
integer i,j,k,ib |
| 1727 |
real taurs,tauoz,tausto,ssatau,asysto,tauto,ssato,asyto |
| 1728 |
real taux,reff1,reff2,g1,g2,asycl(m,n,np) |
| 1729 |
real td(m,n,np+1,2),rr(m,n,np+1,2),tt(m,n,np+1,2), |
| 1730 |
* rs(m,n,np+1,2),ts(m,n,np+1,2) |
| 1731 |
real upflux(m,n,np+1),dwflux(m,n,np+1), |
| 1732 |
* rssab(m,n,np+1),rabx(m,n,np+1),rsabx(m,n,np+1) |
| 1733 |
real fall(m,n,np+1),fclr(m,n,np+1),fsdir(m,n),fsdif(m,n) |
| 1734 |
real asyclt(m,n) |
| 1735 |
real rr1t(m,n),tt1t(m,n),td1t(m,n),rs1t(m,n),ts1t(m,n) |
| 1736 |
real rr2t(m,n),tt2t(m,n),td2t(m,n),rs2t(m,n),ts2t(m,n) |
| 1737 |
|
| 1738 |
c-----hk is the fractional extra-terrestrial solar flux. |
| 1739 |
c the sum of hk is 0.47074. |
| 1740 |
|
| 1741 |
data hk/.00057, .00367, .00083, .00417, |
| 1742 |
* .00600, .00556, .05913, .39081/ |
| 1743 |
|
| 1744 |
c-----xk is the ozone absorption coefficient. unit: /(cm-atm)stp |
| 1745 |
|
| 1746 |
data xk /30.47, 187.2, 301.9, 42.83, |
| 1747 |
* 7.09, 1.25, 0.0345, 0.0539/ |
| 1748 |
|
| 1749 |
c-----ry is the extinction coefficient for Rayleigh scattering. |
| 1750 |
c unit: /mb. |
| 1751 |
|
| 1752 |
data ry /.00604, .00170, .00222, .00132, |
| 1753 |
* .00107, .00091, .00055, .00012/ |
| 1754 |
|
| 1755 |
c-----aerosol single-scattering albedo and asymmetry factor |
| 1756 |
|
| 1757 |
data ssaal,asyal/nband*0.999,nband*0.850/ |
| 1758 |
|
| 1759 |
c-----coefficients for computing the asymmetry factor of ice clouds |
| 1760 |
c from asycl=aig(*,1)+aig(*,2)*reff+aig(*,3)*reff**2 |
| 1761 |
|
| 1762 |
data aig/.74625000,.00105410,-.00000264/ |
| 1763 |
|
| 1764 |
c-----coefficients for computing the asymmetry factor of liquid |
| 1765 |
c clouds from asycl=awg(*,1)+awg(*,2)*reff+awg(*,3)*reff**2 |
| 1766 |
|
| 1767 |
data awg/.82562000,.00529000,-.00014866/ |
| 1768 |
|
| 1769 |
c-----initialize surface reflectances and transmittances |
| 1770 |
|
| 1771 |
do j= 1, n |
| 1772 |
do i= 1, m |
| 1773 |
rr(i,j,np+1,1)=rsuvbm(i,j) |
| 1774 |
rr(i,j,np+1,2)=rsuvbm(i,j) |
| 1775 |
rs(i,j,np+1,1)=rsuvdf(i,j) |
| 1776 |
rs(i,j,np+1,2)=rsuvdf(i,j) |
| 1777 |
td(i,j,np+1,1)=0.0 |
| 1778 |
td(i,j,np+1,2)=0.0 |
| 1779 |
tt(i,j,np+1,1)=0.0 |
| 1780 |
tt(i,j,np+1,2)=0.0 |
| 1781 |
ts(i,j,np+1,1)=0.0 |
| 1782 |
ts(i,j,np+1,2)=0.0 |
| 1783 |
enddo |
| 1784 |
enddo |
| 1785 |
|
| 1786 |
c-----compute cloud asymmetry factor for a mixture of |
| 1787 |
c liquid and ice particles. unit of reff is micrometers. |
| 1788 |
|
| 1789 |
do k= 1, np |
| 1790 |
|
| 1791 |
do j= 1, n |
| 1792 |
do i= 1, m |
| 1793 |
|
| 1794 |
asyclt(i,j)=1.0 |
| 1795 |
|
| 1796 |
taux=taucld(i,j,k,1)+taucld(i,j,k,2) |
| 1797 |
if (taux.gt.0.05 .and. fcld(i,j,k).gt.0.01) then |
| 1798 |
|
| 1799 |
reff1=min(reff(i,j,k,1),130.) |
| 1800 |
reff2=min(reff(i,j,k,2),20.0) |
| 1801 |
|
| 1802 |
g1=(aig(1)+(aig(2)+aig(3)*reff1)*reff1)*taucld(i,j,k,1) |
| 1803 |
g2=(awg(1)+(awg(2)+awg(3)*reff2)*reff2)*taucld(i,j,k,2) |
| 1804 |
asyclt(i,j)=(g1+g2)/taux |
| 1805 |
|
| 1806 |
endif |
| 1807 |
|
| 1808 |
enddo |
| 1809 |
enddo |
| 1810 |
|
| 1811 |
do j=1,n |
| 1812 |
do i=1,m |
| 1813 |
asycl(i,j,k)=asyclt(i,j) |
| 1814 |
enddo |
| 1815 |
enddo |
| 1816 |
|
| 1817 |
enddo |
| 1818 |
|
| 1819 |
do j=1,n |
| 1820 |
do i=1,m |
| 1821 |
fdiruv(i,j)=0.0 |
| 1822 |
fdifuv(i,j)=0.0 |
| 1823 |
enddo |
| 1824 |
enddo |
| 1825 |
|
| 1826 |
c-----integration over spectral bands |
| 1827 |
|
| 1828 |
do 100 ib=1,nband |
| 1829 |
|
| 1830 |
do 300 k= 1, np |
| 1831 |
|
| 1832 |
do j= 1, n |
| 1833 |
do i= 1, m |
| 1834 |
|
| 1835 |
c-----compute ozone and rayleigh optical thicknesses |
| 1836 |
|
| 1837 |
taurs=ry(ib)*dp(i,j,k) |
| 1838 |
tauoz=xk(ib)*oh(i,j,k) |
| 1839 |
|
| 1840 |
c-----compute total optical thickness, single scattering albedo, |
| 1841 |
c and asymmetry factor |
| 1842 |
|
| 1843 |
tausto=taurs+tauoz+taual(i,j,k)+1.0e-8 |
| 1844 |
ssatau=ssaal(ib)*taual(i,j,k)+taurs |
| 1845 |
asysto=asyal(ib)*ssaal(ib)*taual(i,j,k) |
| 1846 |
|
| 1847 |
c-----compute reflectance and transmittance for cloudless layers |
| 1848 |
|
| 1849 |
tauto=tausto |
| 1850 |
ssato=ssatau/tauto+1.0e-8 |
| 1851 |
ssato=min(ssato,0.999999) |
| 1852 |
asyto=asysto/(ssato*tauto) |
| 1853 |
|
| 1854 |
call deledd(tauto,ssato,asyto,csm(i,j), |
| 1855 |
* rr1t(i,j),tt1t(i,j),td1t(i,j)) |
| 1856 |
|
| 1857 |
call sagpol (tauto,ssato,asyto,rs1t(i,j),ts1t(i,j)) |
| 1858 |
|
| 1859 |
c-----compute reflectance and transmittance for cloud layers |
| 1860 |
|
| 1861 |
if (tauclb(i,j,k) .lt. 0.01) then |
| 1862 |
|
| 1863 |
rr2t(i,j)=rr1t(i,j) |
| 1864 |
tt2t(i,j)=tt1t(i,j) |
| 1865 |
td2t(i,j)=td1t(i,j) |
| 1866 |
rs2t(i,j)=rs1t(i,j) |
| 1867 |
ts2t(i,j)=ts1t(i,j) |
| 1868 |
|
| 1869 |
else |
| 1870 |
|
| 1871 |
tauto=tausto+tauclb(i,j,k) |
| 1872 |
ssato=(ssatau+tauclb(i,j,k))/tauto+1.0e-8 |
| 1873 |
ssato=min(ssato,0.999999) |
| 1874 |
asyto=(asysto+asycl(i,j,k)*tauclb(i,j,k))/(ssato*tauto) |
| 1875 |
|
| 1876 |
call deledd(tauto,ssato,asyto,csm(i,j), |
| 1877 |
* rr2t(i,j),tt2t(i,j),td2t(i,j)) |
| 1878 |
|
| 1879 |
tauto=tausto+tauclf(i,j,k) |
| 1880 |
ssato=(ssatau+tauclf(i,j,k))/tauto+1.0e-8 |
| 1881 |
ssato=min(ssato,0.999999) |
| 1882 |
asyto=(asysto+asycl(i,j,k)*tauclf(i,j,k))/(ssato*tauto) |
| 1883 |
|
| 1884 |
call sagpol (tauto,ssato,asyto,rs2t(i,j),ts2t(i,j)) |
| 1885 |
|
| 1886 |
endif |
| 1887 |
|
| 1888 |
enddo |
| 1889 |
enddo |
| 1890 |
|
| 1891 |
do j=1,n |
| 1892 |
do i=1,m |
| 1893 |
rr(i,j,k,1)=rr1t(i,j) |
| 1894 |
enddo |
| 1895 |
enddo |
| 1896 |
do j=1,n |
| 1897 |
do i=1,m |
| 1898 |
tt(i,j,k,1)=tt1t(i,j) |
| 1899 |
enddo |
| 1900 |
enddo |
| 1901 |
do j=1,n |
| 1902 |
do i=1,m |
| 1903 |
td(i,j,k,1)=td1t(i,j) |
| 1904 |
enddo |
| 1905 |
enddo |
| 1906 |
do j=1,n |
| 1907 |
do i=1,m |
| 1908 |
rs(i,j,k,1)=rs1t(i,j) |
| 1909 |
enddo |
| 1910 |
enddo |
| 1911 |
do j=1,n |
| 1912 |
do i=1,m |
| 1913 |
ts(i,j,k,1)=ts1t(i,j) |
| 1914 |
enddo |
| 1915 |
enddo |
| 1916 |
|
| 1917 |
do j=1,n |
| 1918 |
do i=1,m |
| 1919 |
rr(i,j,k,2)=rr2t(i,j) |
| 1920 |
enddo |
| 1921 |
enddo |
| 1922 |
do j=1,n |
| 1923 |
do i=1,m |
| 1924 |
tt(i,j,k,2)=tt2t(i,j) |
| 1925 |
enddo |
| 1926 |
enddo |
| 1927 |
do j=1,n |
| 1928 |
do i=1,m |
| 1929 |
td(i,j,k,2)=td2t(i,j) |
| 1930 |
enddo |
| 1931 |
enddo |
| 1932 |
do j=1,n |
| 1933 |
do i=1,m |
| 1934 |
rs(i,j,k,2)=rs2t(i,j) |
| 1935 |
enddo |
| 1936 |
enddo |
| 1937 |
do j=1,n |
| 1938 |
do i=1,m |
| 1939 |
ts(i,j,k,2)=ts2t(i,j) |
| 1940 |
enddo |
| 1941 |
enddo |
| 1942 |
|
| 1943 |
300 continue |
| 1944 |
|
| 1945 |
c-----flux calculations |
| 1946 |
|
| 1947 |
call cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts, |
| 1948 |
* fclr,fall,fsdir,fsdif) |
| 1949 |
|
| 1950 |
do k= 1, np+1 |
| 1951 |
do j= 1, n |
| 1952 |
do i= 1, m |
| 1953 |
flx(i,j,k)=flx(i,j,k)+fall(i,j,k)*hk(ib) |
| 1954 |
enddo |
| 1955 |
enddo |
| 1956 |
do j= 1, n |
| 1957 |
do i= 1, m |
| 1958 |
flc(i,j,k)=flc(i,j,k)+fclr(i,j,k)*hk(ib) |
| 1959 |
enddo |
| 1960 |
enddo |
| 1961 |
enddo |
| 1962 |
|
| 1963 |
if(ib.eq.nband) then |
| 1964 |
do j=1,n |
| 1965 |
do i=1,m |
| 1966 |
fdirpar(i,j)=fsdir(i,j)*hk(ib) |
| 1967 |
fdifpar(i,j)=fsdif(i,j)*hk(ib) |
| 1968 |
enddo |
| 1969 |
enddo |
| 1970 |
else |
| 1971 |
do j=1,n |
| 1972 |
do i=1,m |
| 1973 |
fdiruv(i,j)=fdiruv(i,j)+fsdir(i,j)*hk(ib) |
| 1974 |
fdifuv(i,j)=fdifuv(i,j)+fsdif(i,j)*hk(ib) |
| 1975 |
enddo |
| 1976 |
enddo |
| 1977 |
endif |
| 1978 |
|
| 1979 |
100 continue |
| 1980 |
|
| 1981 |
return |
| 1982 |
end |
| 1983 |
|
| 1984 |
c********************************************************************* |
| 1985 |
|
| 1986 |
subroutine deledd(tau,ssc,g0,csm,rr,tt,td) |
| 1987 |
|
| 1988 |
c********************************************************************* |
| 1989 |
c |
| 1990 |
c-----uses the delta-eddington approximation to compute the |
| 1991 |
c bulk scattering properties of a single layer |
| 1992 |
c coded following King and Harshvardhan (JAS, 1986) |
| 1993 |
c |
| 1994 |
c inputs: |
| 1995 |
c |
| 1996 |
c tau: the effective optical thickness |
| 1997 |
c ssc: the effective single scattering albedo |
| 1998 |
c g0: the effective asymmetry factor |
| 1999 |
c csm: the effective secant of the zenith angle |
| 2000 |
c |
| 2001 |
c outputs: |
| 2002 |
c |
| 2003 |
c rr: the layer reflection of the direct beam |
| 2004 |
c tt: the layer diffuse transmission of the direct beam |
| 2005 |
c td: the layer direct transmission of the direct beam |
| 2006 |
|
| 2007 |
c********************************************************************* |
| 2008 |
|
| 2009 |
implicit none |
| 2010 |
|
| 2011 |
c-----Explicit Inline Directives |
| 2012 |
|
| 2013 |
#if CRAY |
| 2014 |
#if f77 |
| 2015 |
cfpp$ expand (expmn) |
| 2016 |
#endif |
| 2017 |
#endif |
| 2018 |
real expmn |
| 2019 |
|
| 2020 |
real zero,one,two,three,four,fourth,seven,tumin |
| 2021 |
parameter (one=1., three=3.) |
| 2022 |
parameter (seven=7., two=2.) |
| 2023 |
parameter (four=4., fourth=.25) |
| 2024 |
parameter (zero=0., tumin=1.e-20) |
| 2025 |
|
| 2026 |
c-----input parameters |
| 2027 |
real tau,ssc,g0,csm |
| 2028 |
|
| 2029 |
c-----output parameters |
| 2030 |
real rr,tt,td |
| 2031 |
|
| 2032 |
c-----temporary parameters |
| 2033 |
|
| 2034 |
real zth,ff,xx,taup,sscp,gp,gm1,gm2,gm3,akk,alf1,alf2, |
| 2035 |
* all,bll,st7,st8,cll,dll,fll,ell,st1,st2,st3,st4 |
| 2036 |
c |
| 2037 |
zth = one / csm |
| 2038 |
|
| 2039 |
c delta-eddington scaling of single scattering albedo, |
| 2040 |
c optical thickness, and asymmetry factor, |
| 2041 |
c K & H eqs(27-29) |
| 2042 |
|
| 2043 |
ff = g0*g0 |
| 2044 |
xx = one-ff*ssc |
| 2045 |
taup= tau*xx |
| 2046 |
sscp= ssc*(one-ff)/xx |
| 2047 |
gp = g0/(one+g0) |
| 2048 |
|
| 2049 |
c extinction of the direct beam transmission |
| 2050 |
|
| 2051 |
td = expmn(-taup*csm) |
| 2052 |
|
| 2053 |
c gamma1, gamma2, gamma3 and gamma4. see table 2 and eq(26) K & H |
| 2054 |
c ssc and gp are the d-s single scattering |
| 2055 |
c albedo and asymmetry factor. |
| 2056 |
|
| 2057 |
xx = three*gp |
| 2058 |
gm1 = (seven - sscp*(four+xx))*fourth |
| 2059 |
gm2 = -(one - sscp*(four-xx))*fourth |
| 2060 |
gm3 = (two - zth*xx )*fourth |
| 2061 |
|
| 2062 |
c akk is k as defined in eq(25) of K & H |
| 2063 |
|
| 2064 |
xx = gm1-gm2 |
| 2065 |
akk = sqrt((gm1+gm2)*xx) |
| 2066 |
|
| 2067 |
c alf1 and alf2 are alpha1 and alpha2 from eqs (23) & (24) of K & H |
| 2068 |
|
| 2069 |
alf1 = gm1 - gm3 * xx |
| 2070 |
alf2 = gm2 + gm3 * xx |
| 2071 |
|
| 2072 |
c all is last term in eq(21) of K & H |
| 2073 |
c bll is last term in eq(22) of K & H |
| 2074 |
|
| 2075 |
xx = akk * two |
| 2076 |
all = (gm3 - alf2 * zth )*xx*td |
| 2077 |
bll = (one - gm3 + alf1*zth)*xx |
| 2078 |
|
| 2079 |
xx = akk * zth |
| 2080 |
st7 = one - xx |
| 2081 |
st8 = one + xx |
| 2082 |
|
| 2083 |
xx = akk * gm3 |
| 2084 |
cll = (alf2 + xx) * st7 |
| 2085 |
dll = (alf2 - xx) * st8 |
| 2086 |
|
| 2087 |
xx = akk * (one-gm3) |
| 2088 |
fll = (alf1 + xx) * st8 |
| 2089 |
ell = (alf1 - xx) * st7 |
| 2090 |
|
| 2091 |
st3 = max(abs(st7*st8),tumin) |
| 2092 |
st3 = sign (st3,st7) |
| 2093 |
|
| 2094 |
st2 = expmn(-akk*taup) |
| 2095 |
st4 = st2 * st2 |
| 2096 |
|
| 2097 |
st1 = sscp / ((akk+gm1 + (akk-gm1)*st4) * st3) |
| 2098 |
|
| 2099 |
c rr is r-hat of eq(21) of K & H |
| 2100 |
c tt is diffuse part of t-hat of eq(22) of K & H |
| 2101 |
|
| 2102 |
rr = ( cll-dll*st4 -all*st2)*st1 |
| 2103 |
tt = - ((fll-ell*st4)*td-bll*st2)*st1 |
| 2104 |
|
| 2105 |
rr = max(rr,zero) |
| 2106 |
tt = max(tt,zero) |
| 2107 |
|
| 2108 |
return |
| 2109 |
end |
| 2110 |
|
| 2111 |
c********************************************************************* |
| 2112 |
|
| 2113 |
subroutine sagpol(tau,ssc,g0,rll,tll) |
| 2114 |
|
| 2115 |
c********************************************************************* |
| 2116 |
c-----transmittance (tll) and reflectance (rll) of diffuse radiation |
| 2117 |
c follows Sagan and Pollock (JGR, 1967). |
| 2118 |
c also, eq.(31) of Lacis and Hansen (JAS, 1974). |
| 2119 |
c |
| 2120 |
c-----input parameters: |
| 2121 |
c |
| 2122 |
c tau: the effective optical thickness |
| 2123 |
c ssc: the effective single scattering albedo |
| 2124 |
c g0: the effective asymmetry factor |
| 2125 |
c |
| 2126 |
c-----output parameters: |
| 2127 |
c |
| 2128 |
c rll: the layer reflection of diffuse radiation |
| 2129 |
c tll: the layer transmission of diffuse radiation |
| 2130 |
c |
| 2131 |
c********************************************************************* |
| 2132 |
|
| 2133 |
implicit none |
| 2134 |
|
| 2135 |
c-----Explicit Inline Directives |
| 2136 |
|
| 2137 |
#if CRAY |
| 2138 |
#if f77 |
| 2139 |
cfpp$ expand (expmn) |
| 2140 |
#endif |
| 2141 |
#endif |
| 2142 |
real expmn |
| 2143 |
|
| 2144 |
real one,three,four |
| 2145 |
parameter (one=1., three=3., four=4.) |
| 2146 |
|
| 2147 |
c-----output parameters: |
| 2148 |
|
| 2149 |
real tau,ssc,g0 |
| 2150 |
|
| 2151 |
c-----output parameters: |
| 2152 |
|
| 2153 |
real rll,tll |
| 2154 |
|
| 2155 |
c-----temporary arrays |
| 2156 |
|
| 2157 |
real xx,uuu,ttt,emt,up1,um1,st1 |
| 2158 |
|
| 2159 |
xx = one-ssc*g0 |
| 2160 |
uuu = sqrt( xx/(one-ssc)) |
| 2161 |
ttt = sqrt( xx*(one-ssc)*three )*tau |
| 2162 |
emt = expmn(-ttt) |
| 2163 |
up1 = uuu + one |
| 2164 |
um1 = uuu - one |
| 2165 |
xx = um1*emt |
| 2166 |
st1 = one / ((up1+xx) * (up1-xx)) |
| 2167 |
rll = up1*um1*(one-emt*emt)*st1 |
| 2168 |
tll = uuu*four*emt *st1 |
| 2169 |
|
| 2170 |
return |
| 2171 |
end |
| 2172 |
|
| 2173 |
c******************************************************************* |
| 2174 |
|
| 2175 |
function expmn(fin) |
| 2176 |
|
| 2177 |
c******************************************************************* |
| 2178 |
c compute exponential for arguments in the range 0> fin > -10. |
| 2179 |
|
| 2180 |
parameter (one=1.0, expmin=-10.0) |
| 2181 |
parameter (e1=1.0, e2=-2.507213e-1) |
| 2182 |
parameter (e3=2.92732e-2, e4=-3.827800e-3) |
| 2183 |
real fin,expmn |
| 2184 |
|
| 2185 |
if (fin .lt. expmin) fin = expmin |
| 2186 |
expmn = ((e4*fin + e3)*fin+e2)*fin+e1 |
| 2187 |
expmn = expmn * expmn |
| 2188 |
expmn = one / (expmn * expmn) |
| 2189 |
|
| 2190 |
return |
| 2191 |
end |
| 2192 |
|
| 2193 |
c******************************************************************* |
| 2194 |
|
| 2195 |
subroutine cldflx (m,n,np,ict,icb,cc,rr,tt,td,rs,ts, |
| 2196 |
* fclr,fall,fsdir,fsdif) |
| 2197 |
|
| 2198 |
c******************************************************************* |
| 2199 |
c compute upward and downward fluxes using a two-stream adding method |
| 2200 |
c following equations (3)-(5) of Chou (1992, JAS). |
| 2201 |
c |
| 2202 |
c clouds are grouped into high, middle, and low clouds which are |
| 2203 |
c assumed randomly overlapped. It involves eight sets of calculations. |
| 2204 |
c In each set of calculations, each atmospheric layer is homogeneous, |
| 2205 |
c either with clouds or without clouds. |
| 2206 |
|
| 2207 |
c input parameters: |
| 2208 |
c m: number of soundings in zonal direction |
| 2209 |
c n: number of soundings in meridional direction |
| 2210 |
c np: number of atmospheric layers |
| 2211 |
c ict: the level separating high and middle clouds |
| 2212 |
c icb: the level separating middle and low clouds |
| 2213 |
c cc: effective cloud covers for high, middle and low clouds |
| 2214 |
c tt: diffuse transmission of a layer illuminated by beam radiation |
| 2215 |
c td: direct beam tranmssion |
| 2216 |
c ts: transmission of a layer illuminated by diffuse radiation |
| 2217 |
c rr: reflection of a layer illuminated by beam radiation |
| 2218 |
c rs: reflection of a layer illuminated by diffuse radiation |
| 2219 |
c |
| 2220 |
c output parameters: |
| 2221 |
c fclr: clear-sky flux (downward minus upward) |
| 2222 |
c fall: all-sky flux (downward minus upward) |
| 2223 |
c fdndir: surface direct downward flux |
| 2224 |
c fdndif: surface diffuse downward flux |
| 2225 |
c*********************************************************************c |
| 2226 |
|
| 2227 |
implicit none |
| 2228 |
|
| 2229 |
c-----input parameters |
| 2230 |
|
| 2231 |
integer m,n,np,ict,icb |
| 2232 |
|
| 2233 |
real rr(m,n,np+1,2),tt(m,n,np+1,2),td(m,n,np+1,2) |
| 2234 |
real rs(m,n,np+1,2),ts(m,n,np+1,2) |
| 2235 |
real cc(m,n,3) |
| 2236 |
|
| 2237 |
c-----temporary array |
| 2238 |
|
| 2239 |
integer i,j,k,ih,im,is |
| 2240 |
real rra(m,n,np+1,2,2),tta(m,n,np+1,2,2),tda(m,n,np+1,2,2) |
| 2241 |
real rsa(m,n,np+1,2,2),rxa(m,n,np+1,2,2) |
| 2242 |
real ch(m,n),cm(m,n),ct(m,n),flxdn(m,n,np+1) |
| 2243 |
real fdndir(m,n),fdndif(m,n),fupdif |
| 2244 |
real denm,xx |
| 2245 |
|
| 2246 |
c-----output parameters |
| 2247 |
|
| 2248 |
real fclr(m,n,np+1),fall(m,n,np+1) |
| 2249 |
real fsdir(m,n),fsdif(m,n) |
| 2250 |
|
| 2251 |
c-----initialize all-sky flux (fall) and surface downward fluxes |
| 2252 |
|
| 2253 |
do k=1,np+1 |
| 2254 |
do j=1,n |
| 2255 |
do i=1,m |
| 2256 |
fall(i,j,k)=0.0 |
| 2257 |
enddo |
| 2258 |
enddo |
| 2259 |
enddo |
| 2260 |
|
| 2261 |
do j=1,n |
| 2262 |
do i=1,m |
| 2263 |
fsdir(i,j)=0.0 |
| 2264 |
fsdif(i,j)=0.0 |
| 2265 |
enddo |
| 2266 |
enddo |
| 2267 |
|
| 2268 |
c-----compute transmittances and reflectances for a composite of |
| 2269 |
c layers. layers are added one at a time, going down from the top. |
| 2270 |
c tda is the composite transmittance illuminated by beam radiation |
| 2271 |
c tta is the composite diffuse transmittance illuminated by |
| 2272 |
c beam radiation |
| 2273 |
c rsa is the composite reflectance illuminated from below |
| 2274 |
c by diffuse radiation |
| 2275 |
c tta and rsa are computed from eqs. (4b) and (3b) of Chou |
| 2276 |
|
| 2277 |
c-----for high clouds. indices 1 and 2 denote clear and cloudy |
| 2278 |
c situations, respectively. |
| 2279 |
|
| 2280 |
do 10 ih=1,2 |
| 2281 |
|
| 2282 |
do j= 1, n |
| 2283 |
do i= 1, m |
| 2284 |
tda(i,j,1,ih,1)=td(i,j,1,ih) |
| 2285 |
tta(i,j,1,ih,1)=tt(i,j,1,ih) |
| 2286 |
rsa(i,j,1,ih,1)=rs(i,j,1,ih) |
| 2287 |
tda(i,j,1,ih,2)=td(i,j,1,ih) |
| 2288 |
tta(i,j,1,ih,2)=tt(i,j,1,ih) |
| 2289 |
rsa(i,j,1,ih,2)=rs(i,j,1,ih) |
| 2290 |
enddo |
| 2291 |
enddo |
| 2292 |
|
| 2293 |
do k= 2, ict-1 |
| 2294 |
do j= 1, n |
| 2295 |
do i= 1, m |
| 2296 |
denm = ts(i,j,k,ih)/( 1.-rsa(i,j,k-1,ih,1)*rs(i,j,k,ih)) |
| 2297 |
tda(i,j,k,ih,1)= tda(i,j,k-1,ih,1)*td(i,j,k,ih) |
| 2298 |
tta(i,j,k,ih,1)= tda(i,j,k-1,ih,1)*tt(i,j,k,ih) |
| 2299 |
* +(tda(i,j,k-1,ih,1)*rr(i,j,k,ih) |
| 2300 |
* *rsa(i,j,k-1,ih,1)+tta(i,j,k-1,ih,1))*denm |
| 2301 |
rsa(i,j,k,ih,1)= rs(i,j,k,ih)+ts(i,j,k,ih) |
| 2302 |
* *rsa(i,j,k-1,ih,1)*denm |
| 2303 |
tda(i,j,k,ih,2)= tda(i,j,k,ih,1) |
| 2304 |
tta(i,j,k,ih,2)= tta(i,j,k,ih,1) |
| 2305 |
rsa(i,j,k,ih,2)= rsa(i,j,k,ih,1) |
| 2306 |
enddo |
| 2307 |
enddo |
| 2308 |
enddo |
| 2309 |
|
| 2310 |
c-----for middle clouds |
| 2311 |
|
| 2312 |
do 10 im=1,2 |
| 2313 |
|
| 2314 |
do k= ict, icb-1 |
| 2315 |
do j= 1, n |
| 2316 |
do i= 1, m |
| 2317 |
denm = ts(i,j,k,im)/( 1.-rsa(i,j,k-1,ih,im)*rs(i,j,k,im)) |
| 2318 |
tda(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*td(i,j,k,im) |
| 2319 |
tta(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*tt(i,j,k,im) |
| 2320 |
* +(tda(i,j,k-1,ih,im)*rr(i,j,k,im) |
| 2321 |
* *rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm |
| 2322 |
rsa(i,j,k,ih,im)= rs(i,j,k,im)+ts(i,j,k,im) |
| 2323 |
* *rsa(i,j,k-1,ih,im)*denm |
| 2324 |
enddo |
| 2325 |
enddo |
| 2326 |
enddo |
| 2327 |
|
| 2328 |
10 continue |
| 2329 |
|
| 2330 |
c-----layers are added one at a time, going up from the surface. |
| 2331 |
c rra is the composite reflectance illuminated by beam radiation |
| 2332 |
c rxa is the composite reflectance illuminated from above |
| 2333 |
c by diffuse radiation |
| 2334 |
c rra and rxa are computed from eqs. (4a) and (3a) of Chou |
| 2335 |
|
| 2336 |
c-----for the low clouds |
| 2337 |
|
| 2338 |
do 20 is=1,2 |
| 2339 |
|
| 2340 |
do j= 1, n |
| 2341 |
do i= 1, m |
| 2342 |
rra(i,j,np+1,1,is)=rr(i,j,np+1,is) |
| 2343 |
rxa(i,j,np+1,1,is)=rs(i,j,np+1,is) |
| 2344 |
rra(i,j,np+1,2,is)=rr(i,j,np+1,is) |
| 2345 |
rxa(i,j,np+1,2,is)=rs(i,j,np+1,is) |
| 2346 |
enddo |
| 2347 |
enddo |
| 2348 |
|
| 2349 |
do k=np,icb,-1 |
| 2350 |
do j= 1, n |
| 2351 |
do i= 1, m |
| 2352 |
denm=ts(i,j,k,is)/( 1.-rs(i,j,k,is)*rxa(i,j,k+1,1,is) ) |
| 2353 |
rra(i,j,k,1,is)=rr(i,j,k,is)+(td(i,j,k,is) |
| 2354 |
* *rra(i,j,k+1,1,is)+tt(i,j,k,is)*rxa(i,j,k+1,1,is))*denm |
| 2355 |
rxa(i,j,k,1,is)= rs(i,j,k,is)+ts(i,j,k,is) |
| 2356 |
* *rxa(i,j,k+1,1,is)*denm |
| 2357 |
rra(i,j,k,2,is)=rra(i,j,k,1,is) |
| 2358 |
rxa(i,j,k,2,is)=rxa(i,j,k,1,is) |
| 2359 |
enddo |
| 2360 |
enddo |
| 2361 |
enddo |
| 2362 |
|
| 2363 |
c-----for middle clouds |
| 2364 |
|
| 2365 |
do 20 im=1,2 |
| 2366 |
|
| 2367 |
do k= icb-1,ict,-1 |
| 2368 |
do j= 1, n |
| 2369 |
do i= 1, m |
| 2370 |
denm=ts(i,j,k,im)/( 1.-rs(i,j,k,im)*rxa(i,j,k+1,im,is) ) |
| 2371 |
rra(i,j,k,im,is)= rr(i,j,k,im)+(td(i,j,k,im) |
| 2372 |
* *rra(i,j,k+1,im,is)+tt(i,j,k,im)*rxa(i,j,k+1,im,is))*denm |
| 2373 |
rxa(i,j,k,im,is)= rs(i,j,k,im)+ts(i,j,k,im) |
| 2374 |
* *rxa(i,j,k+1,im,is)*denm |
| 2375 |
enddo |
| 2376 |
enddo |
| 2377 |
enddo |
| 2378 |
|
| 2379 |
20 continue |
| 2380 |
|
| 2381 |
c-----integration over eight sky situations. |
| 2382 |
c ih, im, is denotes high, middle and low cloud groups. |
| 2383 |
|
| 2384 |
do 100 ih=1,2 |
| 2385 |
|
| 2386 |
c-----clear portion |
| 2387 |
|
| 2388 |
if(ih.eq.1) then |
| 2389 |
do j=1,n |
| 2390 |
do i=1,m |
| 2391 |
ch(i,j)=1.0-cc(i,j,1) |
| 2392 |
enddo |
| 2393 |
enddo |
| 2394 |
|
| 2395 |
else |
| 2396 |
|
| 2397 |
c-----cloudy portion |
| 2398 |
|
| 2399 |
do j=1,n |
| 2400 |
do i=1,m |
| 2401 |
ch(i,j)=cc(i,j,1) |
| 2402 |
enddo |
| 2403 |
enddo |
| 2404 |
|
| 2405 |
endif |
| 2406 |
|
| 2407 |
do 100 im=1,2 |
| 2408 |
|
| 2409 |
c-----clear portion |
| 2410 |
|
| 2411 |
if(im.eq.1) then |
| 2412 |
|
| 2413 |
do j=1,n |
| 2414 |
do i=1,m |
| 2415 |
cm(i,j)=ch(i,j)*(1.0-cc(i,j,2)) |
| 2416 |
enddo |
| 2417 |
enddo |
| 2418 |
|
| 2419 |
else |
| 2420 |
|
| 2421 |
c-----cloudy portion |
| 2422 |
|
| 2423 |
do j=1,n |
| 2424 |
do i=1,m |
| 2425 |
cm(i,j)=ch(i,j)*cc(i,j,2) |
| 2426 |
enddo |
| 2427 |
enddo |
| 2428 |
|
| 2429 |
endif |
| 2430 |
|
| 2431 |
do 100 is=1,2 |
| 2432 |
|
| 2433 |
c-----clear portion |
| 2434 |
|
| 2435 |
if(is.eq.1) then |
| 2436 |
|
| 2437 |
do j=1,n |
| 2438 |
do i=1,m |
| 2439 |
ct(i,j)=cm(i,j)*(1.0-cc(i,j,3)) |
| 2440 |
enddo |
| 2441 |
enddo |
| 2442 |
|
| 2443 |
else |
| 2444 |
|
| 2445 |
c-----cloudy portion |
| 2446 |
|
| 2447 |
do j=1,n |
| 2448 |
do i=1,m |
| 2449 |
ct(i,j)=cm(i,j)*cc(i,j,3) |
| 2450 |
enddo |
| 2451 |
enddo |
| 2452 |
|
| 2453 |
endif |
| 2454 |
|
| 2455 |
c-----add one layer at a time, going down. |
| 2456 |
|
| 2457 |
do k= icb, np |
| 2458 |
do j= 1, n |
| 2459 |
do i= 1, m |
| 2460 |
denm = ts(i,j,k,is)/( 1.-rsa(i,j,k-1,ih,im)*rs(i,j,k,is) ) |
| 2461 |
tda(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*td(i,j,k,is) |
| 2462 |
tta(i,j,k,ih,im)= tda(i,j,k-1,ih,im)*tt(i,j,k,is) |
| 2463 |
* +(tda(i,j,k-1,ih,im)*rr(i,j,k,is) |
| 2464 |
* *rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm |
| 2465 |
rsa(i,j,k,ih,im)= rs(i,j,k,is)+ts(i,j,k,is) |
| 2466 |
* *rsa(i,j,k-1,ih,im)*denm |
| 2467 |
enddo |
| 2468 |
enddo |
| 2469 |
enddo |
| 2470 |
|
| 2471 |
c-----add one layer at a time, going up. |
| 2472 |
|
| 2473 |
do k= ict-1,1,-1 |
| 2474 |
do j= 1, n |
| 2475 |
do i= 1, m |
| 2476 |
denm =ts(i,j,k,ih)/(1.-rs(i,j,k,ih)*rxa(i,j,k+1,im,is)) |
| 2477 |
rra(i,j,k,im,is)= rr(i,j,k,ih)+(td(i,j,k,ih) |
| 2478 |
* *rra(i,j,k+1,im,is)+tt(i,j,k,ih)*rxa(i,j,k+1,im,is))*denm |
| 2479 |
rxa(i,j,k,im,is)= rs(i,j,k,ih)+ts(i,j,k,ih) |
| 2480 |
* *rxa(i,j,k+1,im,is)*denm |
| 2481 |
enddo |
| 2482 |
enddo |
| 2483 |
enddo |
| 2484 |
|
| 2485 |
c-----compute fluxes following eq (5) of Chou (1992) |
| 2486 |
|
| 2487 |
c fdndir is the direct downward flux |
| 2488 |
c fdndif is the diffuse downward flux |
| 2489 |
c fupdif is the diffuse upward flux |
| 2490 |
|
| 2491 |
do k=2,np+1 |
| 2492 |
do j=1, n |
| 2493 |
do i=1, m |
| 2494 |
denm= 1./(1.- rxa(i,j,k,im,is)*rsa(i,j,k-1,ih,im)) |
| 2495 |
fdndir(i,j)= tda(i,j,k-1,ih,im) |
| 2496 |
xx = tda(i,j,k-1,ih,im)*rra(i,j,k,im,is) |
| 2497 |
fdndif(i,j)= (xx*rsa(i,j,k-1,ih,im)+tta(i,j,k-1,ih,im))*denm |
| 2498 |
fupdif= (xx+tta(i,j,k-1,ih,im)*rxa(i,j,k,im,is))*denm |
| 2499 |
flxdn(i,j,k)=fdndir(i,j)+fdndif(i,j)-fupdif |
| 2500 |
enddo |
| 2501 |
enddo |
| 2502 |
enddo |
| 2503 |
|
| 2504 |
do j=1, n |
| 2505 |
do i=1, m |
| 2506 |
flxdn(i,j,1)=1.0-rra(i,j,1,im,is) |
| 2507 |
enddo |
| 2508 |
enddo |
| 2509 |
|
| 2510 |
c-----summation of fluxes over all (eight) sky situations. |
| 2511 |
|
| 2512 |
do k=1,np+1 |
| 2513 |
do j=1,n |
| 2514 |
do i=1,m |
| 2515 |
if(ih.eq.1 .and. im.eq.1 .and. is.eq.1) then |
| 2516 |
fclr(i,j,k)=flxdn(i,j,k) |
| 2517 |
endif |
| 2518 |
fall(i,j,k)=fall(i,j,k)+flxdn(i,j,k)*ct(i,j) |
| 2519 |
enddo |
| 2520 |
enddo |
| 2521 |
enddo |
| 2522 |
|
| 2523 |
do j=1,n |
| 2524 |
do i=1,m |
| 2525 |
fsdir(i,j)=fsdir(i,j)+fdndir(i,j)*ct(i,j) |
| 2526 |
fsdif(i,j)=fsdif(i,j)+fdndif(i,j)*ct(i,j) |
| 2527 |
enddo |
| 2528 |
enddo |
| 2529 |
|
| 2530 |
100 continue |
| 2531 |
|
| 2532 |
return |
| 2533 |
end |
| 2534 |
|
| 2535 |
c***************************************************************** |
| 2536 |
|
| 2537 |
subroutine flxco2(m,n,np,swc,swh,csm,df) |
| 2538 |
|
| 2539 |
c***************************************************************** |
| 2540 |
|
| 2541 |
c-----compute the reduction of clear-sky downward solar flux |
| 2542 |
c due to co2 absorption. |
| 2543 |
|
| 2544 |
implicit none |
| 2545 |
|
| 2546 |
c-----input parameters |
| 2547 |
|
| 2548 |
integer m,n,np |
| 2549 |
real csm(m,n),swc(m,n,np+1),swh(m,n,np+1),cah(22,19) |
| 2550 |
|
| 2551 |
c-----output (undated) parameter |
| 2552 |
|
| 2553 |
real df(m,n,np+1) |
| 2554 |
|
| 2555 |
c-----temporary array |
| 2556 |
|
| 2557 |
integer i,j,k,ic,iw |
| 2558 |
real xx,clog,wlog,dc,dw,x1,x2,y2 |
| 2559 |
|
| 2560 |
c******************************************************************** |
| 2561 |
c-----include co2 look-up table |
| 2562 |
|
| 2563 |
include 'cah.dat' |
| 2564 |
save cah |
| 2565 |
|
| 2566 |
c******************************************************************** |
| 2567 |
c-----table look-up for the reduction of clear-sky solar |
| 2568 |
c radiation due to co2. The fraction 0.0343 is the |
| 2569 |
c extraterrestrial solar flux in the co2 bands. |
| 2570 |
|
| 2571 |
do k= 2, np+1 |
| 2572 |
do j= 1, n |
| 2573 |
do i= 1, m |
| 2574 |
xx=1./.3 |
| 2575 |
clog=log10(swc(i,j,k)*csm(i,j)) |
| 2576 |
wlog=log10(swh(i,j,k)*csm(i,j)) |
| 2577 |
ic=int( (clog+3.15)*xx+1.) |
| 2578 |
iw=int( (wlog+4.15)*xx+1.) |
| 2579 |
if(ic.lt.2)ic=2 |
| 2580 |
if(iw.lt.2)iw=2 |
| 2581 |
if(ic.gt.22)ic=22 |
| 2582 |
if(iw.gt.19)iw=19 |
| 2583 |
dc=clog-float(ic-2)*.3+3. |
| 2584 |
dw=wlog-float(iw-2)*.3+4. |
| 2585 |
x1=cah(1,iw-1)+(cah(1,iw)-cah(1,iw-1))*xx*dw |
| 2586 |
x2=cah(ic-1,iw-1)+(cah(ic-1,iw)-cah(ic-1,iw-1))*xx*dw |
| 2587 |
y2=x2+(cah(ic,iw-1)-cah(ic-1,iw-1))*xx*dc |
| 2588 |
if (x1.lt.y2) x1=y2 |
| 2589 |
df(i,j,k)=df(i,j,k)+0.0343*(x1-y2) |
| 2590 |
enddo |
| 2591 |
enddo |
| 2592 |
enddo |
| 2593 |
|
| 2594 |
return |
| 2595 |
end |
Parent Directory
| 
Revision Log
| 
Revision Graph