1!------------------------------------------------------------------------------- 2 3! This file is part of Code_Saturne, a general-purpose CFD tool. 4! 5! Copyright (C) 1998-2021 EDF S.A. 6! 7! This program is free software; you can redistribute it and/or modify it under 8! the terms of the GNU General Public License as published by the Free Software 9! Foundation; either version 2 of the License, or (at your option) any later 10! version. 11! 12! This program is distributed in the hope that it will be useful, but WITHOUT 13! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14! FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 15! details. 16! 17! You should have received a copy of the GNU General Public License along with 18! this program; if not, write to the Free Software Foundation, Inc., 51 Franklin 19! Street, Fifth Floor, Boston, MA 02110-1301, USA. 20 21!------------------------------------------------------------------------------- 22!> \file rayso.f90 23!> 24!> \brief Compute solar fluxes for both clear and cloudy atmosphere following 25!> Lacis and Hansen (1974). The multiple diffusion is taken into account by an 26!> addition method and overlapping between water vapor and liquid water with k 27!> distribution method. 28!> Some improvements from original version concerns: 29!> - introduction of cloud fraction with hazardous recovering 30!> - introduction of aerosol diffusion in the same way as for cloud droplets 31!> but with specific optical properties for aerosols. 32!------------------------------------------------------------------------------- 33! Arguments 34!______________________________________________________________________________. 35! mode name role 36!______________________________________________________________________________! 37!> \param[in] ivertc index of vertical profile 38!> \param[in] k1 index for ground level 39!> \param[in] kmray vertical levels number for radiation 40!> \param[in] heuray Universal time (Hour) 41!> \param[in] imer1 sea index 42!> \param[in] albe albedo 43!> \param[in] qqv optical depth for water vapor (0,z) 44!> \param[in] qqqv idem for intermediate levels 45!> \param[in] qqvinf idem qqv but for altitude above 11000m 46!> \param[in] zqq vertical levels 47!> \param[in] zray altitude (physical mesh) 48!> \param[in] qvray specific umidity for water vapor 49!> \param[in] qlray specific humidity for liquid water 50!> \param[in] fneray cloud fraction 51!> \param[in] romray air density 52!> \param[in] preray pressure 53!> \param[in] temray temperature 54!> \param[in] aeroso aerosol concentration in micro-g/m3 55!> \param[out] fos global downward solar flux at the ground 56!> \param[out] rayst flux divergence of solar radiation 57!_______________________________________________________________________________ 58 59subroutine rayso & 60 (ivertc, k1, kmray, heuray, imer1, albe, & 61 qqv, qqqv, qqvinf, zqq, & 62 zray, qvray, qlray, fneray, & 63 romray, preray, temray, aeroso, fos, rayst, ncray) 64 65!=============================================================================== 66! Module files 67!=============================================================================== 68 69use paramx 70use numvar 71use optcal 72use cstphy 73use entsor 74use parall 75use period 76use ppppar 77use ppthch 78use ppincl 79use cs_c_bindings 80use mesh 81use field 82use atincl, only: kmx, nbmett, sigc, squant, xlat, xlon, sold, & 83 solu, piaero_o3,piaero_h2o, & 84 black_carbon_frac,zaero, gaero_o3, gaero_h2o, & 85 aod_h2o_tot, aod_o3_tot 86use ctincl, only: cp_a, cp_v 87use cstnum, only: epzero, pi 88use radiat 89 90!=============================================================================== 91 92implicit none 93 94! Arguments 95 96integer ivertc, k1, kmray, imer1 97double precision albe, heuray, fos 98double precision qqv(kmx+1), qqqv(kmx+1), qqvinf, zqq(kmx+1) 99double precision qlray(kmx), fneray(kmx), zray(kmx) 100double precision qvray(kmx), preray(kmx) 101double precision aeroso(kmx)!TODO remove 102double precision rayst(kmx), romray(kmx) 103double precision temray(kmx) 104double precision ncray(kmx) 105 106! Local variables 107 108integer i,k,n,l,inua,k1p1,iaer,iaero_top 109integer itop,ibase,itopp1,itopp2,ibasem1 110integer ifac, iz1, iz2, f_id, c_id, iel 111double precision muzero,fo,rr1,m,mbar,rabar,rabar2,rbar 112double precision rabarc,rbarc, refx, trax, refx0, trax0 113double precision qqvtot,y,ystar 114double precision zqm1,zq,xm1,x,xstar,xstarm1 115double precision rrbar,rrbar2s 116double precision tauctot,wh2ol,rm,req,deltaz 117double precision pioc,zbas,dud1 118double precision gasym,drt,tln,drtt1,dtrb1 119double precision kn(8),pkn(8),dqqv 120double precision cphum,qureel 121double precision taua(kmx+1),tauatot 122double precision rabara,rbara 123double precision tauca(kmx+1,8) 124double precision gama1,gama2,kt,gas,fas 125double precision omega, var, zent 126double precision cpvcpa 127! For postprecessing 128double precision soil_direct_flux , soil_global_flux 129double precision soil_direct_flux_h2o, soil_global_flux_h2o 130double precision soil_direct_flux_o3, soil_global_flux_o3 131 132double precision, allocatable:: fabsh2o(:),fabso3(:),tauc(:) 133double precision, allocatable:: tau(:,:),pic(:,:),ref(:,:) 134double precision, allocatable:: reft(:,:),trat(:,:),refts(:,:) 135double precision, allocatable:: refs(:,:),tras(:,:),trats(:,:) 136double precision, allocatable:: refb(:,:),trab(:,:),upw(:,:) 137double precision, allocatable:: refbs(:,:),fabso3c(:,:),tra(:,:) 138double precision, allocatable:: dow(:,:),atln(:,:),absn(:,:) 139double precision, allocatable :: ckup(:), ckdown_r(:), ckdown_f(:) 140double precision, allocatable:: fnebmax(:),fneba(:) 141 142double precision, allocatable, dimension(:,:) :: dowd, trad, trard, ddfso3c 143double precision, allocatable, dimension(:) :: dffsh2o, dffso3 144double precision, allocatable, dimension(:) :: ddfsh2o, ddfso3 145double precision, allocatable, dimension(:) :: drfs, dffs, ddfs, dddsh2o, dddso3 146 147double precision, allocatable, dimension(:) :: dfsh2o, ufsh2o 148double precision, allocatable, dimension(:) :: dfso3, ufso3 149double precision, allocatable, dimension(:,:) :: dfso3c, ufso3c 150double precision, allocatable, dimension(:) :: dfs, ufs 151double precision, dimension(:,:), pointer :: bpro_rad_inc 152double precision, dimension(:,:), pointer :: cpro_ck_up 153double precision, dimension(:,:), pointer :: cpro_ck_down 154! For computing albedo PIC, PIOC 155double precision epsc 156double precision pioco3,pioch2o,gasymo3,gasymh2o 157double precision pic_o3(kmx+1),pic_h2o(kmx+1),gco3(kmx+1),gch2o(kmx+1) 158double precision pioco3_1, pioco3_2,pioch2o_1, pioch2o_2 159double precision rayst_h2o(kmx),rayst_o3(kmx) 160double precision tabara, tabarc 161! 5 minor gas (NO, NO2, CO2, O2, CO) 162double precision Tmg,amg(5),bmg(5),cmg(5),dmg(5),umg(5) 163! For black carbon, 12 bands 164double precision piocv(12), copioc(12), omega0(12), beta1(12) 165double precision beta2(12), beta3(12),beta4(12),copioc20(12) 166double precision nu0, dm0, dm, coeff_E_o3(12), coeff_E_h2o(12) 167double precision pioco3C, pioch2oC 168double precision tauao3(kmx+1) , tauah2o(kmx+1) 169double precision corp, rov 170 171! data for pkn and kn distribution 172data kn/4.d-5,0.002,0.035,0.377,1.95,9.40,44.6,190./ 173data pkn/0.6470,0.0698,0.1443,0.0584,0.0335,0.0225,0.0158,0.0087/ 174 175! Data from Chuang 2002 for calculation of cloud SSA taking into account black carbon 176data beta1/0.2382803,0.2400113,0.2471480,0.2489583,0.2542476,0.2588392,0.2659081,0.2700860,0.2783093,0.2814346,0.2822860,0.1797007/ 177data beta2/0.2940957,0.2936845,0.2880274,0.2871209,0.2824498,0.2775943,0.2698008,0.265296,0.2564840,0.2535739,0.2487382,0.1464709/ 178data beta3/61.83657,58.25082,52.79042,50.06907,45.75322,42.43440,37.03823,32.32349,25.99426,20.05043,12.76966,3.843661/ 179data beta4/574.2111,565.0809,519.0711,494.0088,448.3519,409.9063,348.9051,297.9909,233.7397,175.4385,112.8208,39.24047/ 180data omega0/5189239.d-11,2261712.d-11,1264190.d-11,9446845.d-11,6090293.d-12,3794524.d-12,1735499.d-12,1136807.d-12,2261422.d-12,& 181 1858815.d-11,5551822.d-9,2325124.d-7/ 182data coeff_E_o3/0.0,0.0,0.0,0.0,0.029193795335,0.045219606416,0.16880411522,0.186215099078,0.5705673839,0.0,0.0,0.0/ 183data coeff_E_h2o/0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.27068650951091927, 0.6844296138881737, 0.044883876600907306/ 184 185! Data for calculation of Tmg - Transmission function for minor gases 186data amg/0.0721d0,0.0062d0,0.0326d0,0.0192d0,0.0003d0/ 187data bmg/377.890d0,243.670d0,107.413d0,166.095d0,476.934d0/ 188data cmg/0.5855d0,0.4246d0,0.5501d0,0.4221d0,0.4892d0/ 189data dmg/3.1709d0,1.7222d0,0.9093d0,0.7186d0,0.1261d0/ 190data umg/390.0d0,0.075d0,0.28d0,1.6d0,209500.d0/ 191 192!======================================================================== 193 194allocate(fabsh2o(kmx+1),fabso3(kmx+1),tauc(kmx+1)) 195allocate(tau(kmx+1,8),pic(kmx+1,8),ref(kmx+1,8)) 196allocate(reft(kmx+1,8),trat(kmx+1,8),refts(kmx+1,8)) 197allocate(refs(kmx+1,8),tras(kmx+1,8),trats(kmx+1,8)) 198allocate(refb(kmx+1,8),trab(kmx+1,8),upw(kmx+1,8)) 199allocate(refbs(kmx+1,8),fabso3c(kmx+1,2),tra(kmx+1,8)) 200allocate(dow(kmx+1,8),atln(kmx+1,8),absn(kmx+1,8)) 201allocate(fnebmax(kmx+1),fneba(kmx+1)) 202 203allocate(dfso3c(kmx+1,2), ufso3c(kmx+1,2), ddfso3c(kmx+1,2)) 204allocate(dowd(kmx+1,8), trad(kmx+1,8), trard(kmx+1,8)) 205 206allocate(dfsh2o(kmx+1), ufsh2o(kmx+1)) 207allocate(dfso3(kmx+1), ufso3(kmx+1)) 208allocate(dfs(kmx+1), ufs(kmx+1)) 209allocate(ckup(kmx), ckdown_r(kmx), ckdown_f(kmx)) 210 211allocate(dffsh2o(kmx+1), dffso3(kmx+1)) 212allocate(ddfsh2o(kmx+1), ddfso3(kmx+1)) 213allocate(drfs(kmx+1), dffs(kmx+1), ddfs(kmx+1), dddsh2o(kmx+1), dddso3(kmx+1)) 214 215! 1 - local initializations 216! =========================== 217 218cpvcpa = cp_v / cp_a 219 220inua = 0! TODO test 1 tout le temps 221iaer = 1 ! has aerosols, always, remove 222ibase = 0 223epsc=1.d-8 224do k = 1,kmray 225 if(qlray(k).gt.epsc) inua = 1 226 227 ! TODO usefull ? 228 if (qlray(k).lt.epsc) then 229 qlray(k) = 0.d0 230 fneray(k)=0.d0 231 endif 232enddo 233 234do k = 1, kmx+1 235 fabsh2o(k) = 0.d0 236 fabso3(k) = 0.d0 237 tauc(k) = 0.d0 238 taua(k) = 0.d0 239 tauao3(k)=0.d0 240 tauah2o(k)=0.d0 241 fnebmax(k) = 0.d0 242 gco3(k)=0.d0 243 gch2o(k)=0.d0 244 pic_o3(k)=0.d0 245 pic_h2o(k)=0.d0 246 if(iaer.ge.1) then 247 fneba(k) = 1.d0 248 endif 249 do l = 1, 2 250 fabso3c(k,l) = 0.d0 251 enddo 252 do l = 1, 8 253 tau(k,l) = 0.d0 254 tauca(k,l)=0.d0 255 pic(k,l) = 0.d0 256 atln(k,l) = 0.d0 257 absn(k,l) = 0.d0 258 ref(k,l) = 0.d0 259 reft(k,l) = 0.d0 260 refts(k,l) = 0.d0 261 refs(k,l) = 0.d0 262 refb(k,l) = 0.d0 263 refbs(k,l) = 0.d0 264 trat(k,l) = 0.d0 265 tras(k,l) = 0.d0 266 trab(k,l) = 0.d0 267 trats(k,l) = 0.d0 268 tra(k,l) = 0.d0 269 upw(k,l) = 0.d0 270 dow(k,l) = 0.d0 271 enddo 272enddo 273 274do k = 1, kmx+1 275 do l = 1, 8 276 trad(k,l) = 0.d0 277 dowd(k,l) = 0.d0 278 enddo 279enddo 280refx=0.d0 281trax=0.d0 282refx0=0.d0 283trax0=0.d0 284!initialisation variables for multiple diffusion 285drt=0.d0 286gas=0.d0 287fas=0.d0 288gama1=0.d0 289gama2=0.d0 290kt=0.d0 291tln=0.d0 292! Leighton 1980 293 294 295! id of the top of the aerosol layer 296iaero_top=0 297k1p1 = k1+1 298 299!initialisation Chuang calculations for black carbon in droplets 300dm0=20.d0 !micrometres 301nu0=1.d-8 302!Initialisation of data used for the calculation of SSA using Chuang 2002 303pioco3C=0.d0 304pioch2oC=0.d0 305do k = 1, 12 306 piocv(k)=0.d0 307 copioc(k)=0.d0 308 copioc20(k)=0.d0 309enddo 310 311! 2 - calculation for muzero and solar constant fo 312! =================================================== 313! muzero = cosin of zenithal angle 314! fo = solar constant in watt/m2 315 316! 317! careful : 0. h < heuray < 24. h 318! --------- 319 320qureel = float(squant) 321call raysze(xlat, xlon, qureel, heuray, imer1, albe, muzero, omega, fo) 322! if muzero is negative, it is night and solar radiation is not 323! computed 324 325if (muzero.gt.epzero) then 326 327 ! Correction for very low zenith angles 328 329 rr1 = 0.1255d-2 330 muzero = rr1/(sqrt(muzero**2 + rr1*(rr1 + 2.d0)) - muzero) 331 332 ! Optical air mass 333 ! cf. Kasten, F., Young, A.T., 1989. Revised optical air mass tables and approximation formula. 334 335 ! Note: old formula 336 ! m = 35.d0/sqrt(1224.d0*muzero*muzero + 1.d0) 337 m = 1.d0/(muzero+0.50572d0*(96.07995d0-180.d0/PI*acos(muzero))**(-1.6364d0)) 338 339 mbar = 1.9d0 340 341 ! 3 - albedos for O3 and Rayleigh diffusion 342 343 rabar = 0.219d0/(1.d0 + 0.816d0*muzero) 344 rabar2 = 0.144d0 345 rbar = rabar + (1.d0 - rabar)*(1.d0 - rabar2)*albe/(1.d0 - rabar2*albe) 346 rrbar = 0.28d0/(1.d0 + 6.43d0*muzero) 347 rrbar2s = 0.0685d0 348 349 ! 4 - addition of one level for solar radiation 350 351 zqq(kmray+1) = 16000.d0 352 qqvtot = qqvinf + qqqv(kmray) 353 qqv(kmray+1) = qqvtot - qqvinf/4.d0 354 355 ! Transmission for minor gases 356 Tmg = 1.d0 357 do i = 1, 5 358 Tmg = Tmg* (1.d0 - (amg(i) * m * umg(i)) / & 359 ((1.d0+bmg(i)*m*umg(i))**cmg(i)+dmg(i)*m*umg(i))) 360 enddo 361 362 ! 5 - Solar radiation calculation for cloudy sky 363 ! In order to take into account cloud fraction, multiple diffusion is achieved 364 ! for both cloudy (index 1) and clear (index 2) sky 365 366 ! 5.1 cloud level determination (top for the top of the higher cloud, 367 ! base for the bottom of the lower cloud) 368 369 itop = 0 370 371 do i = kmray, k1p1, -1 372 if(qlray(i).gt.epsc) then 373 if(itop.eq.0) then 374 itop = i 375 ibase = i 376 else 377 ibase = i 378 endif 379 endif 380 enddo 381 382 ! if itop = 0, there is no cloud but, nevertheless, it is possible to execute 383 ! the adding method for the water vapor (SIR band) only 384 385 if(itop.eq.0) then 386 itop = k1 387 ibase = k1 388 endif 389 itopp1 = itop +1 390 itopp2 = itop +2 391 ibasem1 = ibase -1 392 393 ! 5.2 calculation for optical parameters of clouds and aerosols 394 ! (single scattering albedo, optical depth, radius, asymmetry factor) 395 396 fnebmax(kmray+1) = 0.d0 397 tauctot = 0.d0 398 tauatot = 0.d0 399 do i = kmray, k1p1, -1 400 if((i.gt.ibasem1).and.(i.lt.itopp1)) then 401 ! liquid water density in g/m3 in the layers 402 wh2ol = 1.d3*(romray(i)*qlray(i)) 403 ! mean droplet radius in µm 404 rm = 30.d0*wh2ol + 2.d0 405 ! the max of the mean radius is fixed at 10 µm in the considered domain 406 ! and at 2 µm above 407 if(i.le.nbmett) then 408 rm = min(10.d0,rm) 409 else 410 rm = min(2.d0, rm) 411 endif 412 413 ! Efficient radius 414 if (ncray(i).gt.epsc.and.qlray(i).gt.epsc) then 415 ! Simplification: 416 ! Note Old formula 417 ! req = 1.d6*( (3.d0*romray(i)*qlray(i)) / & 418 ! (4.*pi*1000.*ncray(i)*1.d6))**(1./3.) & 419 ! *dexp(sigc**2) 420 req = 1.d3*( (3.d0*romray(i)*qlray(i)) / & 421 (4.d0*pi*ncray(i)))**(1.d0/3.d0) & 422 *dexp(sigc**2) 423 else 424 req = 1.5d0 * rm 425 endif 426 427 ! Clippling: Climatological limits for effective radius 428 if (req .gt. 20.d0) req=20.d0 429 if (req .lt. 1.d0) req=1.d0 430 431 deltaz = zqq(i+1) - zqq(i) 432 ! Cloud optical thickness 433 if (i.eq.k1p1) deltaz = zqq(i+1) - zray(k1) 434 ! req has to be in µm 435 tauc(i) = 1.5d0 * wh2ol * deltaz / req 436 tauctot = tauctot + tauc(i) 437 else 438 tauc(i) = 0.d0 439 req = 0.d0 440 endif 441 442 ! Calculation of aerosol optical depth AOD 443 fneba(i) = 0.d0 444 if((iaer.eq.1).and.(zqq(i).le.zaero)) then 445 iaero_top=max(i,iaero_top) 446 fneba(i) = 1.d0 447 deltaz = zqq(i+1) - zqq(i) 448 if(i.eq.k1p1) deltaz=zqq(i+1) - zray(k1) 449 ! Distribution of AOD on the vertical 450 ! Note, we used a formula based on concentration before v6.2 451 tauao3(i) = aod_o3_tot*deltaz/zqq(iaero_top+1) 452 tauah2o(i) = aod_h2o_tot*deltaz/zqq(iaero_top+1) 453 endif 454 ! Estimation of the law of the cloud fraction 455 fnebmax(i) = max(fnebmax(i+1),fneray(i)) 456 if(black_carbon_frac .gt. epsc) then 457 ! Calculation of SSA for clouds taking into account black carbon fraction - Chuang 2002 458 dm = req*4.d0/3.d0 !mean diameter 459 do k = 5, 12 !5 to 12 because we there is no energy in the 460 !4 first spectral band defined by Chuang 2002 461 copioc20(k)=omega0(k) & 462 + beta1(k)*(1.d0-dexp(-beta3(k)*(black_carbon_frac-nu0))) & 463 + beta2(k)*(1.d0-dexp(-beta4(k)*(black_carbon_frac-nu0))) 464 copioc(k) = (copioc20(k)*dm/dm0) & 465 / (1.d0 + 1.8d0*copioc20(k)*(dm/dm0 - 1.d0)) 466 piocv(k) = 1.d0 - copioc(k) 467 pioco3C = pioco3C+coeff_E_o3(k)*piocv(k) 468 pioch2oC = pioch2oC+coeff_E_h2o(k)*piocv(k) 469 enddo 470 471 pic_h2o(i)=pioch2oC 472 pic_o3(i)=pioco3C 473 pioco3C=0.d0 474 pioch2oC=0.d0 475 else 476 477 ! Calculation of SSA and Asymmetry factor for clouds using- Nielsen 2014 478 ! Note only the first two bands are taken, the third only is 479 ! approximately 0 480 pioco3_1 =( 1.d0 - 33.d-9*req) 481 pioco3_2 = ( 1.d0 - 1.d-7*req ) 482 pioco3=pioco3_1*0.24d0+pioco3_2*0.76d0 483 484 pioch2o_1 =( 0.99999d0 -149.d-7*req ) 485 pioch2o_2 = ( 0.9985d0 -92.d-5*req ) 486 pioch2o=0.60d0*pioch2o_1+0.40d0*pioch2o_2 487 488 gasymo3 =( 0.868d0 + 14.d-5*req & 489 - 61.d-4*dexp(-0.25*req))*pioco3_1*0.24d0 & 490 + ( 0.868d0 + 25.d-5*req & 491 - 63.d-4*dexp(-0.25*req))*pioco3_2*0.76d0 492 493 gasymh2o = ( 0.867d0 + 31.d-5*req & 494 - 78.d-4*dexp(-0.195d0*req))*0.60d0*pioch2o_1 & 495 + ( 0.864d0 + 54.d-5*req & 496 - 0.133d0*dexp(-0.194d0*req))*0.40d0*pioch2o_2 497 498 gco3(i)=gasymo3 499 gch2o(i)=gasymh2o 500 pic_o3(i)=pioco3 501 pic_h2o(i)=pioch2o 502 endif 503 enddo 504 505 fnebmax(k1) = fnebmax(k1p1) 506 507 tauc(kmray+1) = 0.d0 508 509 ! 5.3 O3 absorption in presence of clouds 510 511 ! Calculation of the different albedos for O3 (LH 74) 512 513 ! Asymmetry factor and SSA for liquid water 514 gasym=gco3(itop) 515 516 pioc=pic_o3(itop) 517 !Calculation for cloudy layers 518 call reftra & 519 (pioc, 0.d0, gasym, 0.d0, tauctot, 0.d0, & 520 refx, trax, epsc, 0.d0) 521 522 rabarc=refx 523 tabarc=trax 524 rbarc = rabarc+tabarc*tabarc*albe/(1.d0 - rabarc*albe) 525 !Calculation for aerosol layers 526 call reftra & 527 (0.d0, piaero_o3, 0.d0, gaero_o3, 0.d0, aod_o3_tot, & 528 refx, trax, epsc,0.d0) 529 rabara=refx 530 tabara=trax 531 532 rbara = rabara+tabara*tabara*albe/(1.d0 - rabara*albe) 533 534 ! in the case there is an aerosol layer above the cloud layer 535 536 if((iaer.eq.1).and.(iaero_top.gt.itop)) then 537 itop=iaero_top 538 itopp1=itop+1 539 rbar=rbara 540 rrbar2s=rabara 541 endif 542 543 ! Calculation above the top of the cloud or the aerosol layer 544 545 do i = itop, kmray !calculation have to start at the first level 546 ! ( itop=1 in the case without aerosols and clouds) 547 zqm1 = zqq(i) 548 549 if (i.eq.k1p1) zqm1 = zray(k1) 550 551 zq = zqq(i+1) 552 xm1 = m*rayuoz(zqm1) 553 x = m*rayuoz(zq) 554 555 ! Calculation of heat and radiation fluxes during cloudy sky 556 zbas = zray(itop) 557 xstarm1 = m*rayuoz(zbas) + mbar*(rayuoz(zbas) - rayuoz(zqm1)) 558 xstar = m*rayuoz(zbas) + mbar*(rayuoz(zbas) - rayuoz(zq)) 559 560 ! Heat 561 fabso3c(i,1) = muzero*fo*(raysoz(xm1) - raysoz(x)) & 562 + rbarc*(raysoz(xstar) - raysoz(xstarm1)) 563 564 ! Direct downward radiation 565 ddfso3c(i,1) = muzero*fo*(0.647d0-rrbar-raysoz(x)) 566 ! Downward radiation 567 dfso3c(i,1) = muzero*fo*(0.647d0-rrbar-raysoz(x)) 568 ! Upward (diffuse) radiation 569 ufso3c(i,1) = muzero*fo*(0.647d0-rrbar-raysoz(xstar))*rabarc 570 571 ! Calculation ofheat and radiation fluxes during Clear sky 572 zbas = zray(k1) 573 xstarm1 = m*rayuoz(zbas) + mbar*(rayuoz(zbas) - rayuoz(zqm1)) 574 xstar = m*rayuoz(zbas) + mbar*(rayuoz(zbas) - rayuoz(zq)) 575 !heat 576 fabso3c(i,2) = muzero*fo*(raysoz(xm1) - raysoz(x)) & 577 +rbar*(raysoz(xstar) - raysoz(xstarm1)) 578 579 580 dfso3c(i,2) = muzero*fo*(0.647d0-rrbar-raysoz(x)) & 581 /(1.d0-rrbar2s*albe) 582 583 ddfso3c(i,2) = muzero*fo*(0.647-rrbar-raysoz(x)) 584 ddfso3c(i,2) = min(ddfso3c(i,2),dfso3c(i,2)) 585 ufso3c(i,2) = muzero*fo*(0.647d0-rrbar-raysoz(xstar))*albe & 586 / (1.d0-rrbar2s*albe) 587 ! Summation depending on cloud fraction 588 fabso3(i) = fnebmax(k1p1)*fabso3c(i,1) + (1.d0-fnebmax(k1p1))*fabso3c(i,2) 589 590 591 dfso3(i) = fnebmax(k1p1)*dfso3c(i,1)+ (1.d0-fnebmax(k1p1))*dfso3c(i,2) 592 ddfso3(i) = fnebmax(k1p1)*ddfso3c(i,1)+ (1.d0-fnebmax(k1p1))*ddfso3c(i,2) 593 ufso3(i) = fnebmax(k1p1)*ufso3c(i,1)+ (1.d0-fnebmax(k1p1))*ufso3c(i,2) 594 595 enddo 596 597 ! Calculation under the top of the cloud or the aerosol layer, the adding 598 ! Method with multiple diffusion is used 599 n = 1 !no O3 overlapping 600 do l=k1p1,itopp1 601 tau(l,n)=tauc(l) + tauao3(l) 602 tauca(l,n)= tauao3(l) 603 gasym=gco3(l) 604 pioc=pic_o3(l) 605 ! In the cloud layers 606 if (tauc(l).gt.epsc) then 607 call reftra & 608 (pioc, piaero_o3, gasym, gaero_o3, tauc(l), tauao3(l), & 609 refx, trax, epsc, 0.d0) 610 611 ref(l,n)=fneray(l)*refx 612 tra(l,n)=fneray(l)*trax 613 614 endif 615 ! In the aerosol layers 616 call reftra & 617 ( 0.d0, piaero_o3, 0.d0, gaero_o3, 0.d0, tauao3(l), & 618 refx, trax, epsc, 0.d0) 619 620 ref(l,n)=ref(l,n)+(1.d0-fneray(l))*refx 621 tra(l,n)=tra(l,n) + (1.d0 -fneray(l))*trax 622 623 refs(l,n)=ref(l,n) 624 tras(l,n)=tra(l,n) 625 626 trard(l,n)=fneray(l)* dexp(-m*(tauc(l)+tauao3(l))) & 627 + (1.d0 - fneray(l))*dexp(-m*tauao3(l)) 628 629 end do 630 ! Top boundary conditions 631 tau(itop+1,n)= 0.d0 632 tra(itop+1,n)= 1.d0 633 trard(itop+1,n)= 1.d0 634 trad(itop+1,n)= trard(itop+1,n) 635 ref(itop+1,n)=0.d0 636 tras(itop+1,n)=tra(itop+1,n) 637 refs(itop+1,n)=ref(itop+1,n) 638 trat(itop+1,n)=tra(itop+1,n) 639 reft(itop+1,n)=ref(itop+1,n) 640 refts(itop+1,n)=refs(itop+1,n) 641 trats(itop+1,n)=tras(itop+1,n) 642 ! Bottom boundary conditions 643 tra(k1,n) = 0.d0 644 trard(k1,n) = 0.d0 645 ref(k1,n) = albe 646 tras(k1,n) = 0.d0 647 refs(k1,n) =0.0 648 649 ! Downward addition of layers 650 do l=itop,k1,-1 651 ! Equations 34 of LH74 652 drtt1=1.d0-refts(l+1,n)*ref(l,n) 653 reft(l,n)=reft(l+1,n)+trat(l+1,n)*ref(l,n) & 654 *trats(l+1,n)/drtt1 655 656 trat(l,n)=trat(l+1,n)*tra(l,n)/drtt1 657 658 ! Trad transmission for direct radiation 659 trad(l,n) = trad(l+1,n)*trard(l,n) 660 if(l.gt.k1) then 661 refts(l,n)=refs(l,n)+tras(l,n)*refts(l+1,n)*tra(l,n)/drtt1 662 trats(l,n)=trats(l+1,n)*tras(l,n)/drtt1 663 end if 664 end do 665 666 ! Upward addition of layers 667 refb(k1,n)=ref(k1,n) 668 refbs(k1,n)=refs(k1,n) 669 do l=k1p1,itop 670 dtrb1=1.d0-refb(l-1,n)*refs(l,n) 671 refb(l,n)=ref(l,n)+tra(l,n)*refb(l-1,n)*tras(l,n)/dtrb1 672 end do 673 674 ! Calculation of upward and downward fluxes and absorption 675 do l=itopp1,k1p1,-1 676 dud1=1.d0-refts(l,n)*refb(l-1,n) 677 if(dud1.gt.1.d-30) then 678 679 upw(l,n)=trat(l,n)*refb(l-1,n)/dud1 680 dow(l,n)=trat(l,n)/dud1 681 else 682 upw(l,n)= trat(l,n)*refb(l-1,n) 683 dow(l,n)= trat(l,n) 684 endif 685 dowd(l,n)= trad(l,n) 686 atln(l,n)=((1.d0-reft(k1,n))+upw(l,n)-dow(l,n)) 687 enddo 688 do l = itop, k1p1, -1 689 absn(l,n) = atln(l,n) - atln(l+1,n) 690 ! Fux divergence 691 !flux coming from the top 692 fabso3(l)=fo*muzero*(0.647d0- raysoz(m*rayuoz(zray(itop))))*absn(l,n) 693 enddo 694 695 do i = k1, itop 696 ! addition of ozone absorption for heating in the layers when adding method is used 697 zqm1 = zqq(i) 698 zq = zqq(i+1) 699 if(i.eq.k1p1) zqm1 = zray(k1) 700 xm1 = m*rayuoz(zqm1) 701 x = m*rayuoz(zq) 702 zbas = zray(k1) 703 xstarm1 = m*rayuoz(zbas) + mbar*(rayuoz(zbas) - rayuoz(zqm1)) 704 xstar = m*rayuoz(zbas) + mbar*(rayuoz(zbas) - rayuoz(zq)) 705 !taking into account ozone absorption in the layers with clouds or aerosols 706 fabso3(i) =fabso3(i)+ muzero*fo*(raysoz(xm1) - raysoz(x) + rbar & 707 *(raysoz(xstar) - raysoz(xstarm1))) 708 ! fluxes calculation taking into account ozone absorption 709 dfso3(i)=muzero*fo*(0.647d0-rrbar-raysoz(x))*dow(i+1,n) 710 ddfso3(i)=muzero*fo*(0.647d0-rrbar-raysoz(x))*dowd(i+1,n) 711 ufso3(i)=muzero*fo*(0.647d0-rrbar-raysoz(xstar))*upw(i+1,n) 712 713 enddo 714 715 ! Calculation of upward flux above cloud or aerosol layers taking into account the upward flux transmitted by cloud or aerosol layers 716 do i = itop+1, kmray 717 zq = zray(i) 718 zbas = zray(itop) 719 xstar = mbar*(rayuoz(zbas) - rayuoz(zq)) 720 ufso3(i) =ufso3(itop)*(1.d0 -raysoz(xstar)) 721 enddo 722 723 ! 6.4 Absorption by water vapor and liquid water 724 725 ! In that case we have to solve multiple diffusion. This is achieved by means 726 ! of the adding method following Lacis et Hansen, 1974 727 ! calculation of reflexivity and tarnsmissivity for each vertical layer 728 729 do n = 1, 8 730 do l = k1p1, kmray 731 gasym=gch2o(l) 732 pioc=pic_h2o(l) 733 dqqv = kn(n)*(qqv(l+1) - qqv(l))/10.d0 734 if (l.eq.k1p1) dqqv = kn(n)*qqv(l+1)/10.d0 735 736 ! In the cloud layers 737 tau(l,n) = tauc(l) + dqqv + tauah2o(l) 738 739 if(qlray(l).ge.epsc) then 740 call reftra & 741 (pioc, piaero_h2o, gasym, gaero_h2o, tauc(l) , tauah2o(l), & 742 refx, trax, epsc, dqqv) 743 744 ref(l,n)=fneray(l)*refx 745 tra(l,n)=fneray(l)*trax + (1.d0 - fneray(l))*dexp(-5.d0*dqqv/3.d0) 746 747 if (iaer.eq.1) then 748 call reftra & 749 (0.d0, piaero_h2o, 0.d0, gaero_h2o, 0.d0 , tauah2o(l), & 750 refx0, trax0, epsc, dqqv) 751 752 ref(l,n) = fneray(l)*refx + (1.d0 - fneray(l))*refx0 753 tra(l,n) = fneray(l)*trax + (1.d0 - fneray(l))*trax0 754 endif 755 756 refs(l,n) = ref(l,n) 757 tras(l,n) = tra(l,n) 758 759 ! trard transmissivity for direct radiation 760 trard(l,n) = fneray(l)*dexp(-m*(dqqv+tauc(l)+tauah2o(l))) & 761 +(1.d0-fneray(l))*dexp(-m*(dqqv+tauah2o(l))) 762 763 else 764 765 ! in the clear sky layers 766 ref(l,n) = 0.d0 767 tra(l,n) = dexp(-5.d0*tau(l,n)/3.d0) 768 refs(l,n) = ref(l,n) 769 tras(l,n) = tra(l,n) 770 771 trard(l,n)=dexp(-m*(dqqv+tauah2o(l))) 772 773 774 if(l.ge.itopp1) tra(l,n) = dexp(-m*tau(l,n)) 775 if (iaer.eq.1) then 776 call reftra & 777 (0.d0, piaero_h2o, 0.d0, gaero_h2o, 0.d0, tauah2o(l), & 778 refx, trax, epsc,dqqv) 779 780 ref(l,n)=fneba(l)*refx 781 tra(l,n)=fneba(l)*trax+(1.d0-fneba(l))*dexp(-5.d0*dqqv/3.d0) 782 783 endif 784 785 endif 786 787 enddo 788 789 tau(kmray+1,n) = 0.d0 790 tra(kmray+1,n) = 1.d0 791 792 ! For direct radiation 793 trard(kmray+1,n) = 1.d0 794 trad(kmray+1,n) = trard(kmray+1,n) 795 796 ref(kmray+1,n) = 0.d0 797 tras(kmray+1,n) = tra(kmray+1,n) 798 refs(kmray+1,n) = ref(kmray+1,n) 799 tra(k1,n) = 0.d0 800 ref(k1,n) = albe 801 tras(k1,n) = 0.d0 802 refs(k1,n) = 0.d0 803 804 805 trat(kmray+1,n) = tra(kmray+1,n) 806 reft(kmray+1,n) = ref(kmray+1,n) 807 refts(kmray+1,n) = refs(kmray+1,n) 808 trats(kmray+1,n) = tras(kmray+1,n) 809 fneray(k1) = 0.d0 810 811 ! For direct radiation 812 trad(kmray+1,n) = trard(kmray+1,n) 813 ! downward addition of layers 814 do l = kmray, k1, -1 815 drtt1 = 1.d0 - refts(l+1,n)*ref(l,n) 816 reft(l,n) = reft(l+1,n) & 817 + trat(l+1,n)*ref(l,n)*trats(l+1,n)/drtt1 818 trat(l,n) = trat(l+1,n)*tra(l,n)/drtt1 819 820 ! trad for direct radiation 821 trad(l,n)=trad(l+1,n)*trard(l,n) 822 823 if(l.gt.k1) then 824 refts(l,n) = refs(l,n) & 825 + tras(l,n)*refts(l+1,n)*tra(l,n)/drtt1 826 trats(l,n) = trats(l+1,n)*tras(l,n)/drtt1 827 endif 828 enddo 829 830 ! upward layer addition 831 refb(k1,n) = ref(k1,n) 832 refbs(k1,n) = refs(k1,n) 833 834 do l = k1p1, kmray 835 dtrb1 = 1.d0 - refb(l-1,n)*refs(l,n) 836 refb(l,n) = ref(l,n) + tra(l,n)*refb(l-1,n)*tras(l,n)/dtrb1 837 enddo 838 839 ! calculation of downward and upward fluxes 840 do l = kmray+1, k1p1, -1 841 dud1 = 1.d0 - refts(l,n)*refb(l-1,n) 842 upw(l,n) = trat(l,n)*refb(l-1,n)/dud1 843 dow(l,n) = trat(l,n)/dud1 844 845 ! downward fluxes for direct radiation 846 dowd(l,n) = trad(l,n) 847 848 !calculation of absorption 849 atln(l,n) = pkn(n)*((1.d0 - reft(k1,n)) + upw(l,n) & 850 - dow(l,n)) 851 enddo 852 853 ! absorption in individual layers 854 do l = kmray, k1p1, -1 855 absn(l,n) = atln(l,n) - atln(l+1,n) 856 enddo 857 enddo 858 859 ! summation over frequencies and estimation of absorption integrated 860 ! on the whole spectrum 861 do l = kmray, k1p1, -1 862 fabsh2o(l) = 0.d0 863 do n = 1, 8 864 fabsh2o(l) = fabsh2o(l)+absn(l,n) 865 enddo 866 fabsh2o(l) = fabsh2o(l)*fo*muzero 867 enddo 868 869 ! 5.5 heating in the layers 870 rayst(k1) = 0.d0 871 rayst_h2o(k1) = 0.d0 872 rayst_o3(k1) = 0.d0 873 do i = k1p1, kmray 874 deltaz = zqq(i+1) - zqq(i) 875 if(i.eq.k1p1) deltaz = zqq(i+1) - zray(k1) 876 cphum = cp0*(1.d0 + (cpvcpa - 1.d0)*qvray(i)) 877 rayst(i) = (fabsh2o(i) + fabso3(i))/deltaz/romray(i)/cphum 878 879 rayst_h2o(i) = (fabsh2o(i))/deltaz/romray(i)/cphum 880 rayst_o3(i) = ( fabso3(i))/deltaz/romray(i)/cphum 881 enddo 882 883 ! 5.6 calculation of solar fluxes 884 ! for global radiation, fd for direct radiation for the water vapor band 885 886 do i = k1, kmray 887 dfsh2o(i) = 0.d0 888 ufsh2o(i) = 0.d0 889 ddfsh2o(i) = 0.d0 890 891 do n = 2,8 892 dfsh2o(i) = dfsh2o(i) + pkn(n)*dow(i+1,n) 893 ufsh2o(i) = ufsh2o(i) + pkn(n)*upw(i+1,n) 894 ddfsh2o(i) = ddfsh2o(i) + pkn(n)*dowd(i+1,n) 895 enddo 896 897 dfsh2o(i) = fo*muzero*dfsh2o(i) 898 ufsh2o(i) = fo*muzero*ufsh2o(i) 899 ddfsh2o(i) = fo*muzero*ddfsh2o(i) 900 enddo 901 902 ! 6. Calculation of solar fluxes For the whole spectrum 903 do k = k1, kmray 904 ! Global (down and up) fluxes 905 dfs(k) = dfsh2o(k) + dfso3(k) 906 ufs(k) = ufsh2o(k) + ufso3(k) 907 908 ! direct radiation mod (sum of vapor water band and O3 band) 909 drfs(k) = ddfsh2o(k)+ddfso3(k) 910 ! diffuse radiation (estmated by difference between global and direct) 911 dddsh2o(k) = dfsh2o(k)-ddfsh2o(k) 912 dddso3(k) = dfso3(k)-ddfso3(k) 913 ddfs(k) = dddsh2o(k)+dddso3(k) 914 915 solu(k,ivertc) = ufs(k) 916 sold(k,ivertc) = dfs(k) 917 enddo 918 ! Mutiplication by transmission function for minor gases 919 do k=k1,kmray 920 dfs(k)=Tmg*dfs(k) 921 drfs(k)=Tmg*drfs(k) 922 ufs(k)=Tmg*ufs(k) 923 enddo 924 ! solar heating of the ground surface by the downward global flux 925 fos=dfs(k1)*(1.d0-albe) 926 927 928 929 soil_direct_flux=drfs(k1) 930 soil_global_flux=dfs(k1) 931 soil_direct_flux_h2o=ddfsh2o(k1) 932 soil_global_flux_h2o=dfsh2o(k1) 933 soil_direct_flux_o3=ddfso3(k1) 934 soil_global_flux_o3=dfso3(k1) 935 936 ! For water vapor without clouds and aerosols downward is only direct, 937 ! upward is only diffuse 938 do k = k1, kmray 939 y = m*(qqvtot - qqv(k)) 940 ystar = m*qqvtot + 5.d0/3.d0*qqv(k) 941 corp = (preray(k) / preray(k1))* preray(k1) / 101300.d0!FIXME /p0? 942 rov = romray(k)*(qvray(k)*corp*sqrt(tkelvi/(temray(k) + tkelvi))) 943 ckup(k) = m*dzyama(ystar,rov)/(0.353d0-raysve(ystar)) 944 ckdown_r(k) =m*dzyama(y,rov)/(0.353d0-raysve(y)) 945 ckdown_f(k) = 0.d0 946 enddo 947 948! if muzero < 0, it is night 949else 950 951 muzero = 0.d0 952 do k = k1, kmray 953 rayst(k) = 0.d0 954 955 rayst_h2o(k) = 0.d0 956 rayst_o3(k) = 0.d0 957 958 solu(k,ivertc) = 0.d0 959 sold(k,ivertc) = 0.d0 960 961 ckup(k) = 0.d0 962 ckdown_r(k) = 0.d0 963 ckdown_f(k) = 0.d0 964 enddo 965 soil_direct_flux = 0.d0 966 soil_global_flux = 0.d0 967 soil_direct_flux_h2o=0.0 968 soil_global_flux_h2o=0.0 969 soil_direct_flux_o3=0.0 970 soil_global_flux_o3=0.0 971endif 972 973! TODO compute it 974 975! Compute Boundary conditions for the 3D (Director diFfuse) Solar radiance 976! at the top of the CFD domain 977! and the absorption coefficients 978call field_get_id_try("spectral_rad_incident_flux", f_id) 979 980if (f_id.ge.0) then 981 call field_get_val_v(f_id, bpro_rad_inc) 982 983 call field_get_val_v_by_name("rad_absorption_coeff_up", cpro_ck_up) 984 call field_get_val_v_by_name("rad_absorption_coeff_down", cpro_ck_down) 985 986 c_id = 0 987 ! Direct Solar (denoted by _r) 988 if (iand(rad_atmo_model, 1).eq.1) then 989 990 c_id = c_id + 1 991 992 ! Store the incident radiation of the 1D model 993 do ifac = 1, nfabor 994 995 ! Interpolate at zent 996 zent = cdgfbo(3, ifac) 997 998 call intprz & 999 (kmray, zqq, & 1000 drfs, zent, iz1, iz2, var ) 1001 1002 ! TODO do not multiply and divide by cos(zenital) = muzero 1003 if (muzero.gt.epzero) then 1004 bpro_rad_inc(c_id, ifac) = pi * var / muzero 1005 else 1006 bpro_rad_inc(c_id, ifac) = 0.d0 1007 endif 1008 enddo 1009 1010 ! Store the (downward) absortion coefficient of the 1D model 1011 do iel = 1, ncel 1012 1013 ! Interpolate at zent 1014 zent = xyzcen(3, iel) 1015 1016 call intprz & 1017 (kmray, zqq, & 1018 ckdown_r, zent, iz1, iz2, var ) 1019 1020 cpro_ck_down(c_id, iel) = var! FIXME factor 3/5 ? 1021 enddo 1022 1023 endif 1024 1025 ! Diffuse solar radiation incident 1026 if (iand(rad_atmo_model, 2).eq.2) then 1027 1028 c_id = c_id + 1 1029 do ifac = 1, nfabor 1030 1031 ! Interpolate at zent 1032 zent = cdgfbo(3, ifac) 1033 1034 call intprz & 1035 (kmray, zqq, & 1036 ddfs, zent, iz1, iz2, var ) 1037 1038 bpro_rad_inc(c_id, ifac) = pi * var 1039 1040 enddo 1041 1042 ! Store the (downward and upward) absortion coefficient of the 1D model 1043 do iel = 1, ncel 1044 1045 ! Interpolate at zent 1046 zent = xyzcen(3, iel) 1047 1048 call intprz & 1049 (kmray, zqq, & 1050 ckdown_f, zent, iz1, iz2, var ) 1051 1052 cpro_ck_down(c_id, iel) = var! FIXME factor 3/5 ? 1053 1054 call intprz & 1055 (kmray, zqq, & 1056 ckup, zent, iz1, iz2, var ) 1057 1058 cpro_ck_up(c_id, iel) = var! FIXME factor 3/5 ? 1059 1060 enddo 1061 1062 endif 1063endif 1064 1065deallocate(fabsh2o,fabso3,tauc) 1066deallocate(tau,pic,ref) 1067deallocate(reft,refts) 1068deallocate(refs,tras,trats) 1069deallocate(refb,trab,upw) 1070deallocate(refbs,fabso3c,tra) 1071deallocate(dow,atln,absn) 1072deallocate(fnebmax,fneba) 1073 1074deallocate(dfso3c, ufso3c, ddfso3c) 1075deallocate(dowd, trad, trard) 1076 1077deallocate(dfsh2o, ufsh2o, dfso3, ufso3, dfs, ufs) 1078deallocate(dffsh2o, dffso3) 1079deallocate(ddfsh2o, ddfso3) 1080deallocate(drfs, dffs, ddfs, dddsh2o, dddso3) 1081deallocate(ckup, ckdown_r, ckdown_f) 1082 1083return 1084 1085!=============================================================================== 1086 1087contains 1088 1089 !----------------------------------------------------------------------------- 1090 1091 !> \brief Computes ozone concentration for a given altitude 1092 1093 !> \param[in] zh altitude 1094 1095 function rayuoz(zh) 1096 1097 !=========================================================================== 1098 1099 implicit none 1100 1101 ! Arguments 1102 1103 double precision, intent(in) :: zh ! absolute altitude 1104 1105 ! Local 1106 1107 double precision :: rayuoz 1108 double precision :: a, b, c 1109 1110 !=========================================================================== 1111 1112 a = 0.4d0 1113 b = 20000.d0 1114 c = 5000.d0 1115 1116 rayuoz = a*(1.d0 + dexp(-b/c))/(1.d0 + dexp((zh-b)/c)) 1117 1118 end function rayuoz 1119 1120 !----------------------------------------------------------------------------- 1121 1122 !> \brief Aborption function of the solar radiation by water vapor 1123 1124 !> \param[in] y optical depth for water vapor 1125 1126 function raysve(y) 1127 1128 !=========================================================================== 1129 1130 implicit none 1131 1132 ! Arguments 1133 1134 double precision, intent(in) :: y ! specific humidity 1135 1136 ! Local 1137 1138 double precision :: raysve 1139 1140 !=========================================================================== 1141 1142 raysve = 0.29d0*y/((1.d0 + 14.15d0*y)**0.635d0 + 0.5925d0*y) 1143 1144 end function raysve 1145 1146 !----------------------------------------------------------------------------- 1147 1148 !> \brief Aborption derivative-function of the solar radiation by water vapor 1149 1150 !> \param[in] y optical depth for water vapor 1151 !> \param[in] dy TODO? 1152 1153 function dzyama(y, dy) 1154 1155 implicit none 1156 1157 ! Arguments 1158 1159 double precision, intent(in):: y, dy 1160 1161 double precision:: num,denum 1162 double precision:: dzyama 1163 1164 num = 14.15d0*0.635d0*dy*(1.0d0 + 14.15d0*y)**(0.635d0-1.0d0) + 0.5925d0*dy 1165 denum = (1.0d0 + 14.15d0*y)**0.635d0 + 0.5925d0*y 1166 dzyama= 0.29d0*dy/denum - 0.29d0*y*num/(denum**2.0d0) 1167 1168 end function dzyama 1169 1170 !----------------------------------------------------------------------------- 1171 1172 !> \brief Aborption function of the solar radiation by ozone 1173 1174 !> \param[in] x optical depth for ozone 1175 1176 function raysoz(x) 1177 1178 !=========================================================================== 1179 1180 implicit none 1181 1182 ! Arguments 1183 1184 double precision, intent(in) :: x 1185 1186 ! Local 1187 1188 double precision :: raysoz 1189 double precision :: ao3vis, ao3uv 1190 1191 !=========================================================================== 1192 1193 ao3vis = 0.02118d0*x/(1.d0 + (0.042d0 + 0.000323d0*x)*x) 1194 ao3uv = 1.082d0*x/(1.d0 + 138.6d0*x)**0.805d0 & 1195 + 0.0658d0*x/(1.d0 + (103.6d0*x)**3) 1196 raysoz = ao3vis + ao3uv 1197 1198 end function raysoz 1199 1200 !----------------------------------------------------------------------------- 1201 1202end subroutine rayso 1203 1204!------------------------------------------------------------------------------- 1205!> \brief Compute reflexion and transmission 1206!------------------------------------------------------------------------------- 1207! Arguments 1208!______________________________________________________________________________. 1209! mode name role 1210!______________________________________________________________________________! 1211!> \param[in] pioc Albedo of simple diffusion for cloud (water) 1212!> \param[in] piaero Albedo of simple diffusion for aerosol 1213!> \param[in] gasym Asymmetry factor for clouds 1214!> \param[in] gaero Asymmetry factor for aerosols 1215!> \param[in] tauc Optical depth for clouds 1216!> \param[in] taua Optical depth for aersols 1217!> \param[out] ref Reflexion 1218!> \param[out] tra Transmission 1219!> \param[in] epsc clipping threshold 1220!> \param[in] dqqv Optical depth for Water vapor 1221!_______________________________________________________________________________ 1222 1223subroutine reftra & 1224 (pioc, piaero, gasym, gaero, tauc, taua, & 1225 ref, tra, epsc,dqqv) 1226 1227 !=========================================================================== 1228 1229 implicit none 1230 1231 ! Arguments 1232 1233 double precision, intent(in) :: pioc, piaero, gasym, gaero 1234 double precision, intent(in) :: tauc, taua, dqqv, epsc 1235 double precision, intent(inout) :: ref, tra 1236 1237 ! Local 1238 double precision ::gas, fas, kt, gama1, gama2, tln 1239 double precision :: drt, extlnp, extlnm 1240 double precision :: pic, tau 1241 1242 !=========================================================================== 1243 1244 tau = tauc +taua + dqqv 1245 ! For 0 optical depth 1246 if (tau .lt. epsc) then 1247 ref = 0.d0 1248 tra = 1.d0 1249 else 1250 1251 ! Pure diffusion atmosphere (pioc=1) 1252 if (pioc.ge.(1.d0-epsc)) then !TODO check .and. (taua .le. epsc)) 1253 gama1=(sqrt(3.d0)/2.d0)*(1.d0-gasym) 1254 ref = gama1*tau/(1.d0+gama1*tau) 1255 tra = 1.d0/(1.d0+gama1*tau) 1256 1257 else 1258 pic =(pioc*tauc+piaero*taua)/tau 1259 ! Pure absorbing atmosphere (pioc=0) 1260 if (pic .lt. epsc) then 1261 gama1=dsqrt(3.d0) 1262 ref = 0.d0 1263 tra = dexp(-gama1*tau) 1264 else 1265 1266 gas=(pioc*tauc*gasym+piaero*taua*gaero)& 1267 /(pic*tau) 1268 1269 fas=gas*gas 1270 tau=(1.d0-pic*fas)*tau 1271 pic=pic*(1.d0-fas)/(1.d0-pic*fas) 1272 gas=(gas-fas)/(1.d0-fas) 1273 gama1=(dsqrt(3.d0)/2.d0)*(2.d0-pic*(1.d0+gas)) 1274 gama2=(dsqrt(3.d0)*pic/2.d0)*(1.d0-gas) 1275 kt=dsqrt(gama1*gama1-gama2*gama2) 1276 tln=kt*tau 1277 extlnp=dexp(tln) 1278 extlnm=dexp(-tln) 1279 drt=(kt+gama1)*extlnp+(kt-gama1)*extlnm 1280 ref=gama2*(extlnp-extlnm)/drt 1281 tra=2.d0*kt/drt 1282 endif 1283 1284 endif 1285 1286 endif 1287 1288end subroutine reftra 1289 1290