1head 4.2; 2access; 3symbols; 4locks; strict; 5comment @# @; 6 7 84.2 9date 2004.07.22.15.07.14; author armnphy; state Exp; 10branches; 11next 4.1; 12 134.1 14date 2004.04.20.19.33.53; author armnphy; state Exp; 15branches; 16next 4.0; 17 184.0 19date 2003.09.19.18.49.18; author armnphy; state Exp; 20branches; 21next 3.9; 22 233.9 24date 2003.06.16.18.49.37; author armnphy; state Exp; 25branches; 26next 3.8; 27 283.8 29date 2003.04.01.21.57.11; author armnbil; state Exp; 30branches; 31next ; 32 33 34desc 35@@ 36 37 384.2 39log 40@bidon 41@ 42text 43@!copyright (C) 2001 MSC-RPN COMM %%%RPNPHY%%% 44*** S/P CLDOPTX4 45* 46#include "phy_macros_f.h" 47 subroutine cldoptx4 (LWC,IWC,neb,T,sig,ps,lat,mg,ml,m,lmx,nk, 48 + pbl,ipbl,dz,sdz,eneb,opdepth,asymg, 49 + tlwp,tiwp,topthw,topthi, 50 + ctp,ctt, 51 + omegav,tauae,istcond,satuco, 52 + cw_rad,ioptix,nostrlwc) 53* 54#include "impnone.cdk" 55* 56 integer lmx,m,nk,istcond,cw_rad,ioptix 57 logical satuco,nostrlwc 58 real ipbl(lmx) 59 real LWC(LMX,nk), IWC(LMX,nk), neb(LMX,nk), t(m,nk), sig(LMX,nk) 60 real ps(LMX),lat(LMX),eneb(LMX,nk),mg(LMX),ml(LMX) 61 real opdepth(LMX,nk),asymg(LMX,nk),omegav(LMX,nk) 62 real tauae(LMX,nk,5),pbl(LMX),dz(LMX,nk),sdz(LMX) 63 real tlwp(lmx),tiwp(lmx),topthw(lmx),topthi(lmx) 64 real ctp(lmx),ctt(lmx) 65* 66*AUTHOR 67* L. Garand and H. Barker (April 1995) 68* 69*REVISION 70* 71* 001 R. Sarrazin and L. Garand (May 1996) - Correct bug for omegav 72* and change tuneopi 73* 002 N. Brunet (Oct 96) Correct bug for mg 74* 003 C. Beaudoin (Jan 98) Eliminate fictitious stratospheric clouds 75* above 50 mb for CONDS condensation option 76* 004 B. Bilodeau and L. Garand (Aug 1999) - Add IWC for interaction 77* with microphysics schemes 78* 005 B. Dugas (April 1999) Never any clouds above 70 mb, but this only 79* when the new input parameter climat is true 80* 006 A. Methot and L. Garand (Jun 2000) - introduce a maximum in the 81* total optical depth 82* 007 A. Methot (May 2000) - modify effective radius relationship 83* with lwc 84* 008 A. Methot and L. Garand (Jun 2000) - introduce Hu and Stamnes 85* parameters for liquid water 86* 009 A. Methot and Mailhot (Jun 2000) - introduce Fu & Liou parameters 87* for ice 88* 010 B. Bilodeau (Mar 2001) - Old cldotpx code as option 89* 011 B. Bilodeau (Nov 2001) - Tuneopi = 0.8 for new optical properties 90* 012 B. Bilodeau (Nov 2002) - Back to old optical properties for ice 91* Lakes treated as land 92* 013 B. Bilodeau, P. Vaillancourt and A. Glazer (Dec 2002) 93* - Calculate ctp and ctt 94* 014 B. Dugas - Rename CLIMAT to NOSTRLWC 95* 96* 015 M. Lepine (March 2003) - CVMG... Replacements 97* 016 D. Talbot (June 2003) - IBM conversion 98* - calls to vsexp routine (from massvp4 library) 99* - calls to exponen4 (to calculate power function '**') 100* - divisions replaced by reciprocals 101* 017 B. Bilodeau (Aug 2003) - exponen4 replaced by vspow 102* 018 B. Bilodeau (Feb 2004) - take ml into account for aerosols 103* when ioptix.lt.2 104* 105* 106*OBJECT 107* computes optical parameters as input to visible and infrared 108* radiation also includes aerosol parameterization 109* Optical parameters refer to entire VIS or IR spectrum 110* but could be extended to several bands matching those 111* of the radiation codes. 112* 113*ARGUMENTS 114* - Output - 115* ENEB cloud amount times emissivity in each layer (0. to 1.) 116* (effective nebulosity to be used by IR code) 117* OPDEPTH layer visible optical depth (dimensionless) 118* ASYMG layer visible asymmetry factor (G in literature, 0. to 1. ) 119* OMEGAV layer visible single scattering albedo (0. to 1.) 120* TAUAE layer aerosol optical depth for VIS code 121* IPBL closest model level matching PBL (LMX) 122* TLWP total integrated liquid water path 123* TIWP total integrated ice water path 124* TOPTHW total integrated optical thickness of water (from TLWP) 125* TOPTHI total integrated optical thickness of ice (from TIWP) 126* 127* -Output - 128* LWC TOTAL (liquid and solid) cloud water content for 129* CONDS and OLDSUND schemes (cw_rad=0). 130* Units : Kg water/Kg air (caution: not in Kg/m3) (LMX,NK) 131* 132* -Input - 133* LWC * TOTAL cloud water content for NEWSUND, CONSUN, EXMO 134* and WARM K-Y condensation schemes (cw_rad=1); 135* * LIQUID water content for MIXED PHASE and 136* COLD K-Y schemes (cw_rad=2). 137* IWC ICE water content in Kg water/Kg air (only if cw_rad=2) 138* 139* -Input - 140* NEB layer cloud amount (0. to 1.) (LMX,NK) 141* T layer temperature (K) (M,NK) 142* SIG sigma levels (0. to 1.) (LMX,NK; local sigma) 143* LAT latitude in radians (-pi/2 to +pi/2) (LMX) 144* PBL height of planetary boundary layer in meters (LMX) 145* DZ work array, on output geometrical thickness (LMX,NK) 146* SDZ work array (LMX) 147* MG ground cover (ocean=0.0,land <= 1.) (LMX) 148* ML fraction of lakes (0.-1.) (LMX) 149* PS surface pressure (N/m2) (LMX) 150* LMX number of profiles to compute 151* M first dimension of temperature (usually LMX) 152* NK number of layers 153* ISTCOND stratiform condensation scheme used; 154* SATUCO .TRUE. if water/ice phase for saturation 155* .FALSE. if water phase only for saturation 156* CW_RAD = 0 if no cloud water content is provided as input; 157* = 1 if total water content is provided (in LWC); 158* = 2 if both liquid and ice water contents are provided 159* separately (in LWC and IWC respectively). 160* CW_RAD is defined in phydebu4, based on ISTCOND. 161* NOSTRLWC .TRUE. removes all liquid water content above 70 mb. 162* This should be removed with an updated IR code 163* IOPTIX parameterizations for cloud optical properties 164* = 1 for simpler condensation schemes 165* = 2 for microphysics schemes 166* 167* 168** 169* 170* 171************************************************************************ 172* AUTOMATIC ARRAYS 173************************************************************************ 174* 175 AUTOMATIC ( CLOUD , REAL , (LMX,NK ) ) 176 AUTOMATIC ( DP , REAL , (LMX,NK ) ) 177 AUTOMATIC ( FRAC , REAL , (LMX,NK ) ) 178 AUTOMATIC ( IWCM , REAL , (LMX,NK ) ) 179 AUTOMATIC ( IWP , REAL , (LMX,NK ) ) 180 AUTOMATIC ( LWCM , REAL , (LMX,NK ) ) 181 AUTOMATIC ( LWP , REAL , (LMX,NK ) ) 182 AUTOMATIC ( TRANSMISSINT , REAL , (LMX,NK ) ) 183 AUTOMATIC ( TOP , LOGICAL , (LMX ) ) 184* 185************************************************************************ 186* 187 integer i,k 188 integer ire(m,nk) 189* logical top 190 real lwcm1, tuneopw 191 real rei, rec_rei, ki, iwcm1, tuneopi, omi, gi, ssai 192 real dp1,dp2,dp3 193 real elsa, emiss 194 real third,rec_grav,rec_180,rec_rgasd 195 real ct,aero,rlat,eps 196 real zz, No 197 real tcel(m,nk),aird(m,nk),rew(m,nk),rec_cdd(m,nk),kw(m,nk) 198 real kwf_ire(m,nk),kvf_ire(m,nk),ssf_ire(m,nk),gwf_ire(m,nk) 199 real ei(m,nk),omw(m,nk),ssaw(m,nk),gw(m,nk) 200 real ew(m,nk) 201 real kwf(3,3), kvf(3,3), ssf(3,3), gwf(3,3) 202 external LIQWC 203* 204#include "consphy.cdk" 205* 206c LWP,IWP are liquid, ice water paths in g/m2 207 data third/0.3333333/ 208c diffusivity factor of Elsasser 209 data elsa/1.66/ 210 data eps/1.e-10/ 211 save third,elsa,eps 212* 213 rec_grav=1./grav 214 rec_rgasd=1./rgasd 215* 216 if (IOPTIX.EQ.2) then 217* 218* ************************************************************ 219* LIQUID WATER parameterization after Hu and Stamnes 220* ----------------------------------------------------------- 221* 222* Parameters for the relationship between the 223* diffusivity factor and equivalent radius 224* at 11um wavelenght ( from table 4) 225* After Hu and Stamnes 1993, JCAM p 728-742 226* 2.5 < radius < 12.5 um 227 kwf(1,1) = -3.88e-05 228 kwf(2,1) = 5.24 229 kwf(3,1) = 140. 230* 12.5 <= radius <= 30 um 231 kwf(1,2) = 794. 232 kwf(2,2) = -.148 233 kwf(3,2) = -423. 234* 30.0 < radius < 60 um 235 kwf(1,3) = 1680. 236 kwf(2,3) = -.966 237 kwf(3,3) = -5.84 238 239* Parameters for the relationship between the 240* optical thickness and equivalent radius 241* at .719 um wavelenght ( from table 1) 242* After Hu and Stamnes 1993, JCAM p 728-742 243* 2.5 < radius < 12.5 um 244 kvf(1,1) = 1810. 245 kvf(2,1) = -1.08 246 kvf(3,1) = 6.85 247* 12.5 <= radius <= 30 um 248 kvf(1,2) = 1700. 249 kvf(2,2) = -1.04 250 kvf(3,2) = 1.04 251* 30.0 < radius < 60 um 252 kvf(1,3) = 978. 253 kvf(2,3) = -.816 254 kvf(3,3) = -9.89 255 256* Parameters for the relationship between the 257* asymmetry factor and equivalent radius 258* at .719 um wavelenght ( from table 2) 259* After Hu and Stamnes 1993, JCAM p 728-742 260* 2.5 < radius < 12.5 um 261 ssf(1,1) = 9.95e-7 262 ssf(2,1) = -.856 263 ssf(3,1) = -4.37e-7 264* 12.5 <= radius <= 30 um 265 ssf(1,2) = 1.88e-7 266 ssf(2,2) = 1.32 267 ssf(3,2) = 3.08e-6 268* 30.0 < radius < 60 um 269 ssf(1,3) = 2.03e-5 270 ssf(2,3) = -.332 271 ssf(3,3) = -4.32e-5 272 273* Parameters for the relationship between the 274* single scatering albedo and equivalent radius 275* at .719 um wavelenght ( from their table 3) 276* After Hu and Stamnes 1993, JCAM p 728-742 277* 2.5 < radius < 12.5 um 278 gwf(1,1) = -.141 279 gwf(2,1) = -.694 280 gwf(3,1) = .889 281* 12.5 <= radius <= 30 um 282 gwf(1,2) = -.157 283 gwf(2,2) = -.782 284 gwf(3,2) = .886 285* 30.0 < radius < 60 um 286 gwf(1,3) = -.214 287 gwf(2,3) = -.916 288 gwf(3,3) = .885 289* 290 endif 291* 292* ************************************************************ 293* MISCELLANEOUS 294* ----------------------------------------------------------- 295* 296* tuning for optical thickness in visible 297* (ONLY VIS IS AFFECTED) 298* this tuning affects outgoing and surface 299* radiation; has little impact on atmosperic 300* heating 301* 302 if (ioptix.eq.1) then ! old optical properties 303 tuneopw = 0.30 304 tuneopi = 0.30 305 else if (ioptix.eq.2) then ! new optical properties (for water only) 306 tuneopw = 0.80 307 tuneopi = 0.30 308 endif 309 310* for maximum of lwc 311* 312 call liqwc(asymg,sig,t,ps,lmx,nk,m,satuco) 313* 314* liquid water content if non available as input 315* 316 if(cw_rad.eq.0) then 317 do 33 k=1,nk 318 do 33 i=1,lmx 319* 320* no clouds allowed above 50mb 321* 322 if (sig(i,k).lt.0.050) then 323 lwc(i,k) = 0. 324 else 325 lwc(i,k)=0.4*asymg(i,k) 326 endif 327 33 continue 328 endif 329 330* never any clouds allowed above 70mb in 331* the "NO STRATOSPHERIC LWC" mode 332* 333 if(nostrlwc)then 334 do 34 k=1,nk 335 do 34 i=1,lmx 336 if (sig(i,k) .lt. 0.070) then 337 lwc(i,k) = 0. 338 endif 339 34 continue 340 endif 341* 342* initialize output fields to zero 343* 344 do i=1,lmx 345 tlwp (i) = 0.0 346 tiwp (i) = 0.0 347 topthw(i) = 0.0 348 topthi(i) = 0.0 349 end do 350* 351* ************************************************************ 352* PRELIMINARY WORK 353* ----------------------------------------------------------- 354* 355 do k=1,nk 356 do I=1,lmx 357* 358 lwcm(i,k) = max(lwc(i,k),0.) 359 360 if (cw_rad.le.1) then 361 iwcm(i,k) = 0.0 362 else 363 iwcm(i,k) = max(iwc(i,k),0.) 364 endif 365* 366 cloud(i,k) = max(neb(i,k),0.) 367* 368 if(istcond.gt.1 .and. istcond.lt.5 ) then 369* 370* the following line is an artificial source of clouds 371* when using the "CONDS" condensation option (harmful 372* in the stratosphere) 373* 374 if ((lwcm(i,k)+iwcm(i,k)) .gt. 1.e-6) then 375 cloud(i,k) = max(cloud(i,k) ,0.01) 376 else 377 cloud(i,k) = 0.0 378 endif 379 endif 380 381 if (cloud(i,k) .lt. 0.01) then 382 lwcm(i,k) = 0. 383 iwcm(i,k) = 0. 384 endif 385* 386* max of cloud 387* 388 cloud(i,k) = min(cloud(i,k),1.) 389* 390 if(cw_rad.gt.0) then 391* 392* normalize water contents to get in-cloud values 393* 394 zz=max(cloud(i,k),0.05) 395 lwcm1=lwcm(i,k)/zz 396 iwcm1=iwcm(i,k)/zz 397* 398* consider diabatic lifting limit 399* when Sundquist schem only 400* 401 if ( istcond.lt.5 ) then 402 lwcm(i,k)=min(lwcm1,asymg(i,k)) 403 iwcm(i,k)=min(iwcm1,asymg(i,k)) 404 else 405 lwcm(i,k)=lwcm1 406 iwcm(i,k)=iwcm1 407 endif 408 endif 409* thickness in sigma 410* 411 dp1=0.5*(sig(i,min(k+1,nk))-sig(i,max(k-1,1))) 412 dp2=0.5*(sig(i,1)+sig(i,2)) 413 dp3=0.5*(1.-sig(i,nk)) 414 if (k .eq. 1) then 415 dp(i,k) = dp2 416 else if (k .eq. nk) then 417 dp(i,k) = dp3 418 else 419 dp(i,k) = dp1 420 endif 421 422 dp(i,k)=max(dp(i,k)*ps(i),0.) 423 424 tcel(i,k)=T(i,k)-TCDK 425 426 end do 427 end do 428* 429* LIQUID vs SOLID WATER PARTITION & 430* LIQUID and SOLID WATER PATHS in g/m2 431* 432* In the following, Frac is the fraction of the 433* cloud/precipitation water in the liquid phase 434* after Rockel et al, Beitr. Atmos. Phys, 1991, p.10 435* 436* When this liquid-solid partition is given by 437* the microphysic schem in used ( cw_rad.eq.2 ), 438* frac=1. 439* 440 if ( cw_rad .lt. 2 ) then 441* 442 CALL VSEXP (frac,-.003102*tcel*tcel,nk*lmx) 443 do k=1,nk 444 do I=1,lmx 445* tcel(i,k)=T(i,k)-TCDK 446 if (tcel(i,k) .ge. 0.) then 447 frac(i,k) = 1.0 448 else 449* frac(i,k) = .0059+.9941*exp(-.003102 * tcel(i,k)*tcel(i,k)) 450 frac(i,k) = .0059+.9941*frac(i,k) 451 endif 452 if (frac(i,k) .lt. 0.01) frac(i,k) = 0. 453 454 IWP(i,k) = (1.-frac(i,k))*lwcm(i,k)*dp(i,k)*rec_grav*1000. 455 LWP(i,k) = frac(i,k)*lwcm(i,k)*dp(i,k)*rec_grav*1000. 456 end do 457 end do 458 else 459 do k=1,nk 460 do I=1,lmx 461 frac(i,k) = 1. 462 IWP(i,k) = iwcm(i,k)*dp(i,k)*rec_grav*1000. 463 LWP(i,k) = frac(i,k)*lwcm(i,k)*dp(i,k)*rec_grav*1000. 464 end do 465 end do 466 endif 467 468* 469* end do 470* end do 471* 472* 473* ***************************************************************** 474* MAIN LOOP FOR MICROPHYSICS SCHEMES ("NEW" OPTICAL PROPERTIES) 475* ----------------------------------------------------------------- 476* 477 if (IOPTIX.EQ.2) then 478* 479 do k=1,nk 480 do I=1,lmx 481* EQUIVALENT SIZE of PARTICLES 482* 483*-------------> determines equivalent size of WATER particles 484* set number of drops per cm**3 100 for water 485* and 500 for land 486 487 if (mg(i) .le. 0.5 .and. ml(i) .le. 0.5) then 488* cdd(i,k) = 100. 489 rec_cdd(i,k) = 1. / 100. 490 else 491* cdd(i,k) = 500. 492 rec_cdd(i,k) = 1. / 500. 493 endif 494* 495* aird is air density in kg/m3 496* 497 aird(i,k) = sig(i,k)*ps(i)/T(i,k)*REC_RGASD 498* 499 end do 500 end do 501* 502 CALL VSPOWN1 (REW,(1.+LWCM*1.e4)*LWCM*frac*aird*rec_cdd, 503 + third, nk*lmx) 504* 505 do k=1,nk 506 do I=1,lmx 507* REW(i,k) = 3000. * ( (1.+LWCM(i,k)*1.e4)* 508* + LWCM(i,k)*frac(i,k)*aird(i,k)*rec_cdd(i,k))**third 509 REW(i,k) = 3000. * REW(i,k) 510 REW(i,k) = max( 2.5,REW(i,k)) 511 REW(i,k) = min( 60., REW(i,k)) 512* 513* determines array index for given REW 514 515 ire(i,k) = 2 516 if ( rew(i,k) .lt. 12.5 ) ire(i,k) = 1 517 if ( rew(i,k) .gt. 30.0 ) ire(i,k) = 3 518* avant d'appeler vspownn il faut definir 519 kwf_ire(i,k) = kwf(2,ire(i,k)) 520 kvf_ire(i,k) = kvf(2,ire(i,k)) 521 ssf_ire(i,k) = ssf(2,ire(i,k)) 522 gwf_ire(i,k) = gwf(2,ire(i,k)) 523* 524 end do 525 end do 526* 527 CALL VSPOWNN (kw,rew,kwf_ire,nk*lmx) 528* 529 do k=1,nk 530 do I=1,lmx 531* 532* water diffusivity after Hu and Stammes, JCAM 1993 533* 534* kw = ( kwf(1,ire(i,k))*( rew(i,k)**kwf(2,ire(i,k)) ) + kwf(3,ire(i,k)) )*1.e-3 535 kw(i,k) = ( kwf(1,ire(i,k)) * kw(i,k) + kwf(3,ire(i,k)) )*1.e-3 536* 537 end do 538 end do 539 540 REI = 15. 541 REC_REI = 1. / 15. 542 KI = .0003 + 1.290 * REC_REI 543 CALL VSEXP (EI,-elsa*ki*IWP,nk*lmx) 544 CALL VSEXP (ENEB,-elsa*ki*IWP-kw*LWP,nk*lmx) 545* 546* on elimine une boucle en faisant ceci plus tot 547 CALL VSPOWNN (omw,rew,kvf_ire,nk*lmx) 548 CALL VSPOWNN (ssaw,rew,ssf_ire,nk*lmx) 549 CALL VSPOWNN (gw,rew,gwf_ire,nk*lmx) 550 SSAI = 1.0 - 1.295E-2 - 1.321E-4 *REI 551* 552* 553 do k=1,nk 554 do I=1,lmx 555* 556* 557* emissivity of ice after Ebert and Curry, JGR-1992,p3833 558* using 11.1 micron parameters as representative 559* of entire spectrum assume equivalent ice particles 560* radius of 25 micron will affect both VIS and 561* IR radiation 562* 563* REI = 15. 564* KI = .0003 + 1.290 / REI 565* EI = 1. -exp(-elsa*ki *IWP(i,k)) 566 EI(i,k) = 1. - EI(i,k) 567* 568* 569* compute combined ice/water cloud emissivity assuming 570* the transmission is the product of the ice and water 571* phase transmissions 572* 573* cloud amount is considered outside of the current 574* main loop due to potential optical depth correction 575* 576* cloud emissivity temporarly in ENEB 577 578* ENEB(i,k) = 1. - exp( - elsa*ki *IWP(i,k) - kw(i,k) * LWP(i,k) ) 579 ENEB(i,k) = 1. - ENEB(i,k) 580 neb(i,k) = cloud(i,k) 581* 582* end do 583* end do 584* 585* 586* do k=1,nk 587* do I=1,lmx 588* OPTICAL THICKNESS 589* 590* water optical thickness: Hu and Stammes, JCAM 1993 591* for .719 um 592* 593* omw =LWP(i,k) * ( kvf(1,ire(i,k))*( rew(i,k)**kvf(2,ire(i,k)) ) + kvf(3,ire(i,k)) )*1.e-3 594 omw(i,k)=LWP(i,k) * ( kvf(1,ire(i,k))*( omw(i,k) ) + kvf(3,ire(i,k)) )*1.e-3 595 omw(i,k)=omw(i,k)*tuneopw 596 597 OMI = min(IWP(i,k)*(3.448E-3+2.431*REC_REI) * tuneopi, 25.) 598 599 OPDEPTH(i,k)= max(omw(i,k) + omi,1.e-10) 600* 601* save integrated quantities for output purposes 602* 603 tlwp (i) = tlwp (i) + lwp(i,k) 604 tiwp (i) = tiwp (i) + iwp(i,k) 605 topthw(i) = topthw(i) + omw(i,k) 606 topthi(i) = topthi(i) + omi 607* 608* 609* SINGLE SCATTERING ALBEDO 610* 611* note that this parameter is very close to one for 612* most of VIS but lowers in near infrared; proper 613* spectral weighting is important because small changes 614* will affect substantially the solar heating outgoing 615* radiance or planetary albedo; It has a smaller influence 616* on the surface flux. A SSA of unity will create division 617* by zero in solar code. 618 619* SSAI = 1.0 - 1.295E-2 - 1.321E-4 *REI 620* 621 622* WATER Single scattering Albedo: Hu and Stammes, JCAM 1993 623* for .719 um 624 625* SSAW = 1.- ( ssf(1,ire(i,k))*( rew(i,k)**ssf(2,ire(i,k)) ) + ssf(3,ire(i,k)) ) 626 SSAW(i,k)= 1.- ( ssf(1,ire(i,k))*( SSAW(i,k) ) + ssf(3,ire(i,k)) ) 627* 628* weighting by optical depth 629* 630 OMEGAV(i,k)= (OMW(i,k) * SSAW(i,k) + OMI* SSAI)/ (OMW(i,k)+OMI+eps) 631 OMEGAV(i,k)=max(OMEGAV(i,k),0.9850) 632 OMEGAV(i,k)=min(OMEGAV(i,k),0.999999) 633* 634* Ice Asymmetry factor: Ebert and Curry JGR 1992 Eq.8 635* for 25 micron particles and a weighting for the 636* five bands given in their Table 2 637* 638 GI = min(0.777 + 5.897E-4 * REI, 0.9999) 639* 640* 641* Water Asymetry factor: Hu and Stammes, JCAM 1993 642* for .719 um 643 644* GW = ( gwf(1,ire(i,k))*( rew(i,k)**gwf(2,ire(i,k)) ) + gwf(3,ire(i,k)) ) 645 GW(i,k) = ( gwf(1,ire(i,k))*( GW(i,k) ) + gwf(3,ire(i,k)) ) 646* 647* weighting by SSA * opt depth 648* 649 ASYMG(i,k) = (SSAI*OMI*GI + SSAW(i,k)*OMW(i,k)*GW(i,k))/(SSAI*OMI+SSAW(i,k)*OMW(i,k)+eps) 650* 651 asymg(i,k) = max(asymg(i,k), 0.75 ) 652* 653 asymg(i,k) = min(asymg(i,k),0.9999) 654* 655* geometrical thickness 656* 657 dz(i,k)= dp(i,k)/aird(i,k)*rec_grav 658* 659* 660 end do 661 end do 662 663 else if (IOPTIX.EQ.1) then 664* 665* **************************************************************** 666* MAIN LOOP FOR CONVENTIONAL CONDENSATION SCHEMES 667* ("OLD" OPTICAL PROPERTIES) 668* ---------------------------------------------------------------- 669* 670 REI = 15. 671 REC_REI = 1. / 15. 672 KI = .0003 + 1.290 * REC_REI 673* 674 CALL VSEXP (EW,-0.087*elsa*LWP,nk*lmx) 675 CALL VSEXP (EI,-elsa*ki*IWP,nk*lmx) 676* 677 do k=1,nk 678 do I=1,lmx 679* 680* emissivity of water after Stephens, JAS 78 681* we take a 0.144 factor as average of his 682* downward and upward emissivity which divided 683* by diffusivity factor elsa yields 0.087 684* 685* EW = 1. -exp(-0.087*elsa* LWP(i,k)) 686 EW(i,k) = 1. - EW(i,k) 687* 688* emissivity of ice after Ebert and Curry, JGR-1992,p3833 689* using 11.1 micron parameters as representative 690* of entire spectrum assume equivalent ice particles 691* radius of 25 micron will affect both VIS and 692* IR radiation 693* 694* EI(i,k) = 1. -exp(-elsa*ki *IWP(i,k)) 695 EI(i,k) = 1. - EI(i,k) 696* 697* compute combined ice/water cloud emissivity assuming 698* the transmission is the product of the ice and water 699* phase transmissions 700* 701 EMISS = 1. - (1.-EI(i,k))* (1.-EW(i,k)) 702* 703* effective cloud cover 704* 705 ENEB(i,k)= cloud(i,k)*emiss 706* 707* black clouds for istcond non 3 708* 709* eneb(i,k)= cvmgt(cloud(i,k),eneb(i,k),istcond.NE.3) 710* 711* optical thickness at nadir is computed 712* set number of drops per cm**3 100 for water 713* and 500 for land 714* 715* The following line is commented out to ensure bitwise 716* validation of the GEMDM global model. It should be 717* activated in a future version of the code. 718 if (mg(i).le.0.5) then 719* cdd(i,k) = 100. 720 rec_cdd(i,k) = 1. / 100. 721 else 722* cdd(i,k) = 500. 723 rec_cdd(i,k) = 1. / 500. 724 endif 725* 726* aird is air density in kg/m3 727 aird(i,k) = sig(i,k)*ps(i)/T(i,k)*REC_RGASD 728* 729 end do 730 end do 731* 732 CALL VSPOWN1 (REW,LWCM*frac*aird*rec_cdd,third,nk*lmx) 733* 734 do k=1,nk 735 do I=1,lmx 736* 737* this parameterization from H. Barker, based 738* on aircraft data range 4-17 micron is that specified 739* by Slingo for parameterizations 740* 741* REW(i,k) = max(4., 754.6 * (LWCM(i,k)*frac(i,k)*aird(i,k)*rec_cdd(i,k))**third) 742 REW(i,k) = max(4., 754.6 * REW(i,k)) 743 REW(i,k) = min(REW(i,k),17.) 744 745* 746* slingo JAS 89 p 1420 for weighted average of 747* bands 1-5 (all spectrum) 748* 749 OMW(i,k) = LWP(i,k)* (2.622E-2 + 1.356/REW(i,k) ) 750* 751* follows Ebert and Curry, 1992 752* no variation as function of band 753* ice optical depth limited to 25; water to 125. 754* 755 omw(i,k)= min(omw(i,k)*tuneopw,125.) 756 OMI = min(IWP(i,k)*(3.448E-3+2.431*REC_REI) * tuneopi, 25.) 757 OPDEPTH(i,k)= max(omw(i,k) + omi,1.e-10) 758* 759* 760* save integrated quantities for output purposes 761* 762 tlwp (i) = tlwp (i) + lwp(i,k) 763 tiwp (i) = tiwp (i) + iwp(i,k) 764 topthw(i) = topthw(i) + omw(i,k) 765 topthi(i) = topthi(i) + omi 766* 767 SSAI = 1.0 - 1.295E-2 - 1.321E-4 *REI 768* Slingo JAS 1989 for weighted average bands 1-4 p 1420 769 SSAW(i,k) = 1.0 - 6.814E-3 - 4.842E-4 *REW(i,k) 770* 771* 772* weighting by optical depth 773* 774 OMEGAV(i,k)= (OMW(i,k) * SSAW(i,k) + OMI* SSAI)/ (OMW(i,k)+OMI+eps) 775 OMEGAV(i,k)=max(OMEGAV(i,k),0.9850) 776 OMEGAV(i,k)=min(OMEGAV(i,k),0.9999) 777* 778* 779* Ice Asymmetry factor: Ebert and Curry JGR 1992 Eq.8 780* for 25 micron particles and a weighting for the 781* five bands given in their Table 2 782* 783 GI = min(0.777 + 5.897E-4 * REI, 0.9999) 784* 785* Water Asymetry factor: Slingo 1989 786* 787 GW(i,k) = min(0.804 + 3.850E-3*REW(i,k), 0.9999) 788* 789* weighting by SSA * opt depth 790* 791 ASYMG(i,k) = (SSAI*OMI*GI + SSAW(i,k)*OMW(i,k)*GW(i,k))/(SSAI*OMI+SSAW(i,k)*OMW(i,k)+eps) 792* 793 asymg(i,k) = max(asymg(i,k),GI) 794* 795 asymg(i,k) = min(asymg(i,k),0.9999) 796* 797* geometrical thickness 798* 799 dz(i,k)= dp(i,k)/aird(i,k)*rec_grav 800* 801 end do 802 end do 803* 804* 805 endif 806* 807* 808* 809* ************************************************************ 810* END OF MAIN LOOPS 811* ----------------------------------------------------------- 812 813* 814 if (IOPTIX.EQ.2) then 815* 816* ****************************************************************** 817* CORRECTION FOR MAXIMUM TOTAL OPTICAL DEPTH & 818* FINAL CLOUD*EMISSIVITY CALCULATIONS 819* ------------------------------------------------------------------ 820* 821* temporary use of ipbl as a work field. 822* ipbl is 1. when there is no correction 823* otherwise ipbl is a number between zero and one 824* 825 do i=1, lmx 826 ipbl(i)=min( 1., 20./( max(topthi(i)+topthw(i), 1.e-10) ) ) 827 enddo 828* only optical depth and cloud emissivity 829* need to be re-scaled 830* 831 do k=1, nk 832 do i=1, lmx 833 OPDEPTH(i,k) = ipbl(i) * OPDEPTH(i,k) 834 eneb(i,k) = neb(i,k) * 835 $ ( 1. - ( 1.-eneb(i,k) )**ipbl(i) ) 836 enddo 837 enddo 838* 839* 840 endif 841* 842* Diagnostics : cloud top pressure (ctp) and temperature (ctt) 843* 844 do i=1,lmx 845 ctp (i) = 110000. 846 ctt (i) = 310. 847 top(i) = .true. 848 transmissint(i,1) = 1. - eneb(i,1) 849 if ( (1.-transmissint(i,1)) .gt. 0.99 .and. top(i) ) then 850 ctp(i) = sig(i,1)*ps(i) 851 ctt(i) = t(i,1) 852 top(i) = .false. 853 end if 854 end do 855 856 do k=2,nk 857 do i=1,lmx 858 transmissint(i,k) = transmissint(i,k-1) * (1.-eneb(i,k)) 859 if ( (1.-transmissint(i,k)) .gt. 0.99 .and. top(i) ) then 860 ctp(i) = sig(i,k)*ps(i) 861 ctt(i) = t(i,k) 862 top(i) = .false. 863 end if 864 end do 865 end do 866* 867* 868* ****************************************************************** 869* AEROSOLS LOADING 870* ------------------------------------------------------------------ 871* 872 do 10 I =1,lmx 873 sdz(i) = 0. 874 ipbl(i) = 0. 875 10 continue 876* 877 do 5 k=nk,1,-1 878 do 6 I=1,lmx 879 if ( int(ipbl(i)).eq.0 ) then 880* pbl heiht recomputed as sum of layer thicknesses 881 sdz(i)=sdz(i)+dz(i,k) 882* level closest to pbl 883 if (sdz(i) .gt. pbl(i)) ipbl(i) = float(k) 884 endif 885 6 continue 886 5 continue 887* 888* distributing aerosols 889* optical thickness higher over land; 890* decreases with latitude 891* See Toon and Pollack, J. Appl. Meteor, 1976, p.235 892* aero being total optical thickness, we distribute it 893* within PBL in proportion with geometrical thickness 894* above pbl aerosol optical depth is kept negligible 895* 896 ct=2./pi 897 REC_180=1./180. 898 do 3 k=1,nk 899 do 4 i=1,lmx 900 tauae(i,k,1)=1.e-10 901 tauae(i,k,2)=1.e-10 902 rlat = lat(i)*pi*REC_180 903 if ( k.ge.INT(ipbl(i)) ) then 904* 905 if ( mg(i).ge.0.5.or.ml(i).ge.0.5) then 906* 907* over land 908* 909 aero = 0.25 - 0.2*ct * abs(lat(i)) 910 tauae(i,k,1)= max(aero * dz(i,k)/sdz(i), 1.e-10) 911 else 912* 913* over ocean 914* 915 aero = 0.13 - 0.1*ct * abs(lat(i)) 916 tauae(i,k,2)= max(aero *dz(i,k)/sdz(i), 1.e-10) 917 endif 918 endif 919* 920* other types of aerosols set to negligible 921* 922 tauae(i,k,3)=1.e-10 923 tauae(i,k,4)=1.e-10 924 tauae(i,k,5)=1.e-10 925 4 continue 926 3 continue 927 return 928 end 929@ 930 931 9324.1 933log 934@bidon 935@ 936text 937@d60 2 938d863 1 939a863 10 940 if ( (ioptix.eq.2 .and. 941 + (mg(i).ge.0.5.or.ml(i).ge.0.5)) .or. 942 + (ioptix.lt.2 .and. 943 + mg(i).ge.0.5 ) ) then 944* 945* The following line is commented out to ensure bitwise 946* validation of the GEMDM global model. It should be 947* activated in a future version of the code, and 948* the previous lines should be deleted. 949c if ( mg(i).ge.0.5.or.ml(i).ge.0.5) then 950@ 951 952 9534.0 954log 955@La version 4.0 de la physique a ete creee le 4 septembre 2003. 956 957C'est la premiere version de la programmatheque des parametrages 958destinee a etre mise en exploitation sur le nouveau calculateur 959IBM. 960 961Principaux changements : 962---------------------- 963 964 * Correction de bogues introduits lors de la conversion vers IBM 965 * Nouvelles optimisations pour IBM 966 * Modifications pour permettre l'interpolation des profils 967 atmospheriques pres de la surface 968@ 969text 970@@ 971 972 9733.9 974log 975@La version 3.9 de la physique a ete creee le 16 juin 2003. 976 977Elle constitue la premiere version de conversion vers le 978calculateur IBM. 979 980Le nouveau code de "gravity wave drag" sgoflx3.ftn est une 981copie du code linearise lin_sgoflx1.ftn. 982@ 983text 984@d59 1 985d458 2 986a459 2 987 CALL EXPONEN4 (REW,(1.+LWCM*1.e4)*LWCM*frac*aird*rec_cdd, 988 + third, nk*lmx, nk*lmx, 1) 989d474 1 990a474 1 991* pour faire fonctionner exponen4 il faut definir 992d483 1 993a483 1 994 CALL EXPONEN4 (kw,rew,kwf_ire,nk*lmx,nk*lmx,nk*lmx) 995d503 3 996a505 3 997 CALL EXPONEN4 (omw,rew,kvf_ire,nk*lmx,nk*lmx,nk*lmx) 998 CALL EXPONEN4 (ssaw,rew,ssf_ire,nk*lmx,nk*lmx,nk*lmx) 999 CALL EXPONEN4 (gw,rew,gwf_ire,nk*lmx,nk*lmx,nk*lmx) 1000a540 4 1001* CALL EXPONEN4 (omw,rew,kvf_ire,nk*lmx,nk*lmx,nk*lmx) 1002* CALL EXPONEN4 (ssaw,rew,ssf_ire,nk*lmx,nk*lmx,nk*lmx) 1003* CALL EXPONEN4 (gw,rew,gwf_ire,nk*lmx,nk*lmx,nk*lmx) 1004* SSAI = 1.0 - 1.295E-2 - 1.321E-4 *REI 1005d688 1 1006a688 1 1007 CALL EXPONEN4 (REW,LWCM*frac*aird*rec_cdd,third,nk*lmx,nk*lmx,1) 1008a786 2 1009 CALL EXPONEN4 (eneb,1.-eneb,ipbl,nk*lmx,nk*lmx,lmx) 1010* 1011d790 2 1012a791 2 1013* eneb(i,k) = neb(i,k) * ( 1. - ( 1.-eneb(i,k) )**ipbl(i) ) 1014 eneb(i,k) = neb(i,k) * ( 1. - eneb(i,k) ) 1015@ 1016 1017 10183.8 1019log 1020@La version 3.8 de la physique a ete creee le 12 mars 2003 1021 1022Principaux changements par rapport a la version 3.72 : 1023---------------------------------------------------- 1024 1025 * contenu de la rustine r2 de la version 3.72 1026 * code developpe pour le modele regional a 15 km 1027 - MOISTKE (refonte) 1028 - MIXPHASE (avec BLG) 1029 - KTRSNT 1030 - proprietes optiques des nuages 1031 * option ADVECTKE reactivee 1032 * BOUJO disponible dans eturbl 1033 * ajouts importants au code de physique linearisee 1034 * nouvelles cles : AS,BETA,KKL,KFCPCP,SHLCVT(2),STRATOS,QCO2 1035 * nombreux diagnostics supplementaires 1036 * optimisation des series temporelles 1037 * diagnostics supplementaires pour CHRONOS et AURAMS 1038 * correction d'une multitude de bogues mineurs 1039@ 1040text 1041@d54 7 1042a124 3 1043#if defined (CVMG) 1044#include "cvmg.cdk" 1045#endif 1046d142 5 1047a146 4 1048 integer i,k, ire 1049 logical top 1050 real rew, kw, lwcm1, tuneopw, omw, gw, ssaw 1051 real rei, ki, iwcm1, tuneopi, omi, gi, ssai 1052d148 2 1053a149 2 1054 real ei, ew, elsa, emiss 1055 real aird, third, tcel, cdd 1056d152 4 1057d168 3 1058a273 1 1059 lwc(i,k)=0.4*asymg(i,k) 1060d277 5 1061a281 1 1062 lwc(i,k)=cvmgt(0.,lwc(i,k),sig(i,k).lt.0.050) 1063d291 3 1064a293 1 1065 lwc(i,k)=cvmgt(0.,lwc(i,k),sig(i,k).lt.0.070) 1066d329 5 1067a333 2 1068 cloud(i,k) = cvmgt(max(cloud(i,k) ,0.01),0.0, 1069 + (lwcm(i,k)+iwcm(i,k)).gt.1.e-6) 1070d336 4 1071a339 2 1072 lwcm(i,k)= cvmgt(0., lwcm(i,k), cloud(i,k) .lt. 0.01) 1073 iwcm(i,k)= cvmgt(0., iwcm(i,k), cloud(i,k) .lt. 0.01) 1074d369 8 1075a376 2 1076 dp(i,k)=cvmgt(dp2,dp1,k.eq.1) 1077 dp(i,k)=cvmgt(dp3,dp(i,k),k.eq.nk) 1078d378 5 1079d397 11 1080a407 4 1081 tcel=T(i,k)-TCDK 1082 frac(i,k) = cvmgt(1.,.0059+.9941*exp(-.003102 * tcel*tcel), 1083 x tcel.ge.0.) 1084 frac(i,k)= cvmgt(0.,frac(i,k), frac(i,k).lt.0.01) 1085d409 4 1086a412 1 1087 IWP(i,k) = (1.-frac(i,k))*lwcm(i,k)*dp(i,k)/grav*1000. 1088d414 2 1089a415 1 1090* 1091d417 4 1092a420 1 1093 IWP(i,k) = iwcm(i,k)*dp(i,k)/grav*1000. 1094a422 1 1095 LWP(i,k) = frac(i,k)*lwcm(i,k)*dp(i,k)/grav*1000. 1096d424 2 1097a425 2 1098 end do 1099 end do 1100d442 7 1101a448 1 1102 cdd =cvmgt(100.,500.,mg(i).le.0.5.and.ml(i).le.0.5) 1103d452 7 1104a458 1 1105 aird = sig(i,k)*ps(i)/T(i,k)/RGASD 1106d460 7 1107a466 4 1108 REW = 3000. * ( (1.+LWCM(i,k)*1.e4)* 1109 + LWCM(i,k)*frac(i,k)*aird/cdd )**third 1110 REW = max( 2.5,REW) 1111 REW = min( 60., REW) 1112d470 13 1113a482 4 1114 ire = 2 1115 if ( rew .lt. 12.5 ) ire = 1 1116 if ( rew .gt. 30.0 ) ire = 3 1117 1118d484 2 1119d488 6 1120d495 16 1121a510 2 1122 kw = ( kwf(1,ire)*( rew**kwf(2,ire) ) + kwf(3,ire) )*1.e-3 1123 1124d518 5 1125a522 3 1126 REI = 15. 1127 KI = .0003 + 1.290 / REI 1128 EI = 1. -exp(-elsa*ki *IWP(i,k)) 1129a523 1 1130 1131d533 2 1132a534 1 1133 ENEB(i,k) = 1. - exp( - elsa*ki *IWP(i,k) - kw * LWP(i,k) ) 1134d537 10 1135d551 4 1136a554 3 1137 1138 omw=LWP(i,k) * ( kvf(1,ire)*( rew**kvf(2,ire) ) + kvf(3,ire) )*1.e-3 1139 omw=omw*tuneopw 1140d556 1 1141a556 1 1142 OMI = min(IWP(i,k)*(3.448E-3+2.431/REI) * tuneopi, 25.) 1143d558 1 1144a558 1 1145 OPDEPTH(i,k)= max(omw + omi,1.e-10) 1146d564 1 1147a564 1 1148 topthw(i) = topthw(i) + omw 1149d578 1 1150a578 1 1151 SSAI = 1.0 - 1.295E-2 - 1.321E-4 *REI 1152d584 2 1153a585 1 1154 SSAW= 1.- ( ssf(1,ire)*( rew**ssf(2,ire) ) + ssf(3,ire) ) 1155d589 1 1156a589 1 1157 OMEGAV(i,k)= (OMW * SSAW + OMI* SSAI)/ (OMW+OMI+eps) 1158d603 2 1159a604 1 1160 GW = ( gwf(1,ire)*( rew**gwf(2,ire) ) + gwf(3,ire) ) 1161d608 1 1162a608 1 1163 ASYMG(i,k) = (SSAI*OMI*GI + SSAW*OMW*GW)/(SSAI*OMI+SSAW*OMW+eps) 1164d616 1 1165a616 1 1166 dz(i,k)= dp(i,k)/(aird*grav) 1167d629 7 1168d644 2 1169a645 1 1170 EW = 1. -exp(-0.087*elsa* LWP(i,k)) 1171d653 2 1172a654 3 1173 REI = 15. 1174 KI = .0003 + 1.290 / REI 1175 EI = 1. -exp(-elsa*ki *IWP(i,k)) 1176d660 1 1177a660 1 1178 EMISS = 1. - (1.-EI)* (1.-EW) 1179d677 7 1180a683 2 1181c cdd =cvmgt(100.,500.,mg(i).le.0.5.and.ml(i).le.0.5) 1182 cdd =cvmgt(100.,500.,mg(i).le.0.5) 1183d686 10 1184a695 2 1185 aird = sig(i,k)*ps(i)/T(i,k)/RGASD 1186 1187d700 3 1188a702 2 1189 REW = max(4., 754.6 * (LWCM(i,k)*frac(i,k)*aird/cdd)**third) 1190 REW = min(REW,17.) 1191d708 1 1192a708 1 1193 OMW = LWP(i,k)* (2.622E-2 + 1.356/REW ) 1194d714 3 1195a716 3 1196 omw= min(omw*tuneopw,125.) 1197 OMI = min(IWP(i,k)*(3.448E-3+2.431/REI) * tuneopi, 25.) 1198 OPDEPTH(i,k)= max(omw + omi,1.e-10) 1199d723 1 1200a723 1 1201 topthw(i) = topthw(i) + omw 1202d728 1 1203a728 1 1204 SSAW = 1.0 - 6.814E-3 - 4.842E-4 *REW 1205d733 1 1206a733 1 1207 OMEGAV(i,k)= (OMW * SSAW + OMI* SSAI)/ (OMW+OMI+eps) 1208d746 1 1209a746 1 1210 GW = min(0.804 + 3.850E-3*REW, 0.9999) 1211d750 1 1212a750 1 1213 ASYMG(i,k) = (SSAI*OMI*GI + SSAW*OMW*GW)/(SSAI*OMI+SSAW*OMW+eps) 1214d758 1 1215a758 1 1216 dz(i,k)= dp(i,k)/(aird*grav) 1217d790 2 1218d795 2 1219a796 2 1220 eneb(i,k) = neb(i,k) * 1221 $ ( 1. - ( 1.-eneb(i,k) )**ipbl(i) ) 1222d844 1 1223a844 1 1224 ipbl(i)=cvmgt(float(k),ipbl(i),sdz(i).gt.pbl(i)) 1225d858 1 1226d863 1 1227a863 1 1228 rlat = lat(i)*pi/180. 1229@ 1230