head 4.2; access; symbols; locks; strict; comment @# @; 4.2 date 2004.07.22.15.07.14; author armnphy; state Exp; branches; next 4.1; 4.1 date 2004.04.20.19.33.53; author armnphy; state Exp; branches; next 4.0; 4.0 date 2003.09.19.18.49.18; author armnphy; state Exp; branches; next 3.9; 3.9 date 2003.06.16.18.49.37; author armnphy; state Exp; branches; next 3.8; 3.8 date 2003.04.01.21.57.11; author armnbil; state Exp; branches; next ; desc @@ 4.2 log @bidon @ text @!copyright (C) 2001 MSC-RPN COMM %%%RPNPHY%%% *** S/P CLDOPTX4 * #include "phy_macros_f.h" subroutine cldoptx4 (LWC,IWC,neb,T,sig,ps,lat,mg,ml,m,lmx,nk, + pbl,ipbl,dz,sdz,eneb,opdepth,asymg, + tlwp,tiwp,topthw,topthi, + ctp,ctt, + omegav,tauae,istcond,satuco, + cw_rad,ioptix,nostrlwc) * #include "impnone.cdk" * integer lmx,m,nk,istcond,cw_rad,ioptix logical satuco,nostrlwc real ipbl(lmx) real LWC(LMX,nk), IWC(LMX,nk), neb(LMX,nk), t(m,nk), sig(LMX,nk) real ps(LMX),lat(LMX),eneb(LMX,nk),mg(LMX),ml(LMX) real opdepth(LMX,nk),asymg(LMX,nk),omegav(LMX,nk) real tauae(LMX,nk,5),pbl(LMX),dz(LMX,nk),sdz(LMX) real tlwp(lmx),tiwp(lmx),topthw(lmx),topthi(lmx) real ctp(lmx),ctt(lmx) * *AUTHOR * L. Garand and H. Barker (April 1995) * *REVISION * * 001 R. Sarrazin and L. Garand (May 1996) - Correct bug for omegav * and change tuneopi * 002 N. Brunet (Oct 96) Correct bug for mg * 003 C. Beaudoin (Jan 98) Eliminate fictitious stratospheric clouds * above 50 mb for CONDS condensation option * 004 B. Bilodeau and L. Garand (Aug 1999) - Add IWC for interaction * with microphysics schemes * 005 B. Dugas (April 1999) Never any clouds above 70 mb, but this only * when the new input parameter climat is true * 006 A. Methot and L. Garand (Jun 2000) - introduce a maximum in the * total optical depth * 007 A. Methot (May 2000) - modify effective radius relationship * with lwc * 008 A. Methot and L. Garand (Jun 2000) - introduce Hu and Stamnes * parameters for liquid water * 009 A. Methot and Mailhot (Jun 2000) - introduce Fu & Liou parameters * for ice * 010 B. Bilodeau (Mar 2001) - Old cldotpx code as option * 011 B. Bilodeau (Nov 2001) - Tuneopi = 0.8 for new optical properties * 012 B. Bilodeau (Nov 2002) - Back to old optical properties for ice * Lakes treated as land * 013 B. Bilodeau, P. Vaillancourt and A. Glazer (Dec 2002) * - Calculate ctp and ctt * 014 B. Dugas - Rename CLIMAT to NOSTRLWC * * 015 M. Lepine (March 2003) - CVMG... Replacements * 016 D. Talbot (June 2003) - IBM conversion * - calls to vsexp routine (from massvp4 library) * - calls to exponen4 (to calculate power function '**') * - divisions replaced by reciprocals * 017 B. Bilodeau (Aug 2003) - exponen4 replaced by vspow * 018 B. Bilodeau (Feb 2004) - take ml into account for aerosols * when ioptix.lt.2 * * *OBJECT * computes optical parameters as input to visible and infrared * radiation also includes aerosol parameterization * Optical parameters refer to entire VIS or IR spectrum * but could be extended to several bands matching those * of the radiation codes. * *ARGUMENTS * - Output - * ENEB cloud amount times emissivity in each layer (0. to 1.) * (effective nebulosity to be used by IR code) * OPDEPTH layer visible optical depth (dimensionless) * ASYMG layer visible asymmetry factor (G in literature, 0. to 1. ) * OMEGAV layer visible single scattering albedo (0. to 1.) * TAUAE layer aerosol optical depth for VIS code * IPBL closest model level matching PBL (LMX) * TLWP total integrated liquid water path * TIWP total integrated ice water path * TOPTHW total integrated optical thickness of water (from TLWP) * TOPTHI total integrated optical thickness of ice (from TIWP) * * -Output - * LWC TOTAL (liquid and solid) cloud water content for * CONDS and OLDSUND schemes (cw_rad=0). * Units : Kg water/Kg air (caution: not in Kg/m3) (LMX,NK) * * -Input - * LWC * TOTAL cloud water content for NEWSUND, CONSUN, EXMO * and WARM K-Y condensation schemes (cw_rad=1); * * LIQUID water content for MIXED PHASE and * COLD K-Y schemes (cw_rad=2). * IWC ICE water content in Kg water/Kg air (only if cw_rad=2) * * -Input - * NEB layer cloud amount (0. to 1.) (LMX,NK) * T layer temperature (K) (M,NK) * SIG sigma levels (0. to 1.) (LMX,NK; local sigma) * LAT latitude in radians (-pi/2 to +pi/2) (LMX) * PBL height of planetary boundary layer in meters (LMX) * DZ work array, on output geometrical thickness (LMX,NK) * SDZ work array (LMX) * MG ground cover (ocean=0.0,land <= 1.) (LMX) * ML fraction of lakes (0.-1.) (LMX) * PS surface pressure (N/m2) (LMX) * LMX number of profiles to compute * M first dimension of temperature (usually LMX) * NK number of layers * ISTCOND stratiform condensation scheme used; * SATUCO .TRUE. if water/ice phase for saturation * .FALSE. if water phase only for saturation * CW_RAD = 0 if no cloud water content is provided as input; * = 1 if total water content is provided (in LWC); * = 2 if both liquid and ice water contents are provided * separately (in LWC and IWC respectively). * CW_RAD is defined in phydebu4, based on ISTCOND. * NOSTRLWC .TRUE. removes all liquid water content above 70 mb. * This should be removed with an updated IR code * IOPTIX parameterizations for cloud optical properties * = 1 for simpler condensation schemes * = 2 for microphysics schemes * * ** * * ************************************************************************ * AUTOMATIC ARRAYS ************************************************************************ * AUTOMATIC ( CLOUD , REAL , (LMX,NK ) ) AUTOMATIC ( DP , REAL , (LMX,NK ) ) AUTOMATIC ( FRAC , REAL , (LMX,NK ) ) AUTOMATIC ( IWCM , REAL , (LMX,NK ) ) AUTOMATIC ( IWP , REAL , (LMX,NK ) ) AUTOMATIC ( LWCM , REAL , (LMX,NK ) ) AUTOMATIC ( LWP , REAL , (LMX,NK ) ) AUTOMATIC ( TRANSMISSINT , REAL , (LMX,NK ) ) AUTOMATIC ( TOP , LOGICAL , (LMX ) ) * ************************************************************************ * integer i,k integer ire(m,nk) * logical top real lwcm1, tuneopw real rei, rec_rei, ki, iwcm1, tuneopi, omi, gi, ssai real dp1,dp2,dp3 real elsa, emiss real third,rec_grav,rec_180,rec_rgasd real ct,aero,rlat,eps real zz, No real tcel(m,nk),aird(m,nk),rew(m,nk),rec_cdd(m,nk),kw(m,nk) real kwf_ire(m,nk),kvf_ire(m,nk),ssf_ire(m,nk),gwf_ire(m,nk) real ei(m,nk),omw(m,nk),ssaw(m,nk),gw(m,nk) real ew(m,nk) real kwf(3,3), kvf(3,3), ssf(3,3), gwf(3,3) external LIQWC * #include "consphy.cdk" * c LWP,IWP are liquid, ice water paths in g/m2 data third/0.3333333/ c diffusivity factor of Elsasser data elsa/1.66/ data eps/1.e-10/ save third,elsa,eps * rec_grav=1./grav rec_rgasd=1./rgasd * if (IOPTIX.EQ.2) then * * ************************************************************ * LIQUID WATER parameterization after Hu and Stamnes * ----------------------------------------------------------- * * Parameters for the relationship between the * diffusivity factor and equivalent radius * at 11um wavelenght ( from table 4) * After Hu and Stamnes 1993, JCAM p 728-742 * 2.5 < radius < 12.5 um kwf(1,1) = -3.88e-05 kwf(2,1) = 5.24 kwf(3,1) = 140. * 12.5 <= radius <= 30 um kwf(1,2) = 794. kwf(2,2) = -.148 kwf(3,2) = -423. * 30.0 < radius < 60 um kwf(1,3) = 1680. kwf(2,3) = -.966 kwf(3,3) = -5.84 * Parameters for the relationship between the * optical thickness and equivalent radius * at .719 um wavelenght ( from table 1) * After Hu and Stamnes 1993, JCAM p 728-742 * 2.5 < radius < 12.5 um kvf(1,1) = 1810. kvf(2,1) = -1.08 kvf(3,1) = 6.85 * 12.5 <= radius <= 30 um kvf(1,2) = 1700. kvf(2,2) = -1.04 kvf(3,2) = 1.04 * 30.0 < radius < 60 um kvf(1,3) = 978. kvf(2,3) = -.816 kvf(3,3) = -9.89 * Parameters for the relationship between the * asymmetry factor and equivalent radius * at .719 um wavelenght ( from table 2) * After Hu and Stamnes 1993, JCAM p 728-742 * 2.5 < radius < 12.5 um ssf(1,1) = 9.95e-7 ssf(2,1) = -.856 ssf(3,1) = -4.37e-7 * 12.5 <= radius <= 30 um ssf(1,2) = 1.88e-7 ssf(2,2) = 1.32 ssf(3,2) = 3.08e-6 * 30.0 < radius < 60 um ssf(1,3) = 2.03e-5 ssf(2,3) = -.332 ssf(3,3) = -4.32e-5 * Parameters for the relationship between the * single scatering albedo and equivalent radius * at .719 um wavelenght ( from their table 3) * After Hu and Stamnes 1993, JCAM p 728-742 * 2.5 < radius < 12.5 um gwf(1,1) = -.141 gwf(2,1) = -.694 gwf(3,1) = .889 * 12.5 <= radius <= 30 um gwf(1,2) = -.157 gwf(2,2) = -.782 gwf(3,2) = .886 * 30.0 < radius < 60 um gwf(1,3) = -.214 gwf(2,3) = -.916 gwf(3,3) = .885 * endif * * ************************************************************ * MISCELLANEOUS * ----------------------------------------------------------- * * tuning for optical thickness in visible * (ONLY VIS IS AFFECTED) * this tuning affects outgoing and surface * radiation; has little impact on atmosperic * heating * if (ioptix.eq.1) then ! old optical properties tuneopw = 0.30 tuneopi = 0.30 else if (ioptix.eq.2) then ! new optical properties (for water only) tuneopw = 0.80 tuneopi = 0.30 endif * for maximum of lwc * call liqwc(asymg,sig,t,ps,lmx,nk,m,satuco) * * liquid water content if non available as input * if(cw_rad.eq.0) then do 33 k=1,nk do 33 i=1,lmx * * no clouds allowed above 50mb * if (sig(i,k).lt.0.050) then lwc(i,k) = 0. else lwc(i,k)=0.4*asymg(i,k) endif 33 continue endif * never any clouds allowed above 70mb in * the "NO STRATOSPHERIC LWC" mode * if(nostrlwc)then do 34 k=1,nk do 34 i=1,lmx if (sig(i,k) .lt. 0.070) then lwc(i,k) = 0. endif 34 continue endif * * initialize output fields to zero * do i=1,lmx tlwp (i) = 0.0 tiwp (i) = 0.0 topthw(i) = 0.0 topthi(i) = 0.0 end do * * ************************************************************ * PRELIMINARY WORK * ----------------------------------------------------------- * do k=1,nk do I=1,lmx * lwcm(i,k) = max(lwc(i,k),0.) if (cw_rad.le.1) then iwcm(i,k) = 0.0 else iwcm(i,k) = max(iwc(i,k),0.) endif * cloud(i,k) = max(neb(i,k),0.) * if(istcond.gt.1 .and. istcond.lt.5 ) then * * the following line is an artificial source of clouds * when using the "CONDS" condensation option (harmful * in the stratosphere) * if ((lwcm(i,k)+iwcm(i,k)) .gt. 1.e-6) then cloud(i,k) = max(cloud(i,k) ,0.01) else cloud(i,k) = 0.0 endif endif if (cloud(i,k) .lt. 0.01) then lwcm(i,k) = 0. iwcm(i,k) = 0. endif * * max of cloud * cloud(i,k) = min(cloud(i,k),1.) * if(cw_rad.gt.0) then * * normalize water contents to get in-cloud values * zz=max(cloud(i,k),0.05) lwcm1=lwcm(i,k)/zz iwcm1=iwcm(i,k)/zz * * consider diabatic lifting limit * when Sundquist schem only * if ( istcond.lt.5 ) then lwcm(i,k)=min(lwcm1,asymg(i,k)) iwcm(i,k)=min(iwcm1,asymg(i,k)) else lwcm(i,k)=lwcm1 iwcm(i,k)=iwcm1 endif endif * thickness in sigma * dp1=0.5*(sig(i,min(k+1,nk))-sig(i,max(k-1,1))) dp2=0.5*(sig(i,1)+sig(i,2)) dp3=0.5*(1.-sig(i,nk)) if (k .eq. 1) then dp(i,k) = dp2 else if (k .eq. nk) then dp(i,k) = dp3 else dp(i,k) = dp1 endif dp(i,k)=max(dp(i,k)*ps(i),0.) tcel(i,k)=T(i,k)-TCDK end do end do * * LIQUID vs SOLID WATER PARTITION & * LIQUID and SOLID WATER PATHS in g/m2 * * In the following, Frac is the fraction of the * cloud/precipitation water in the liquid phase * after Rockel et al, Beitr. Atmos. Phys, 1991, p.10 * * When this liquid-solid partition is given by * the microphysic schem in used ( cw_rad.eq.2 ), * frac=1. * if ( cw_rad .lt. 2 ) then * CALL VSEXP (frac,-.003102*tcel*tcel,nk*lmx) do k=1,nk do I=1,lmx * tcel(i,k)=T(i,k)-TCDK if (tcel(i,k) .ge. 0.) then frac(i,k) = 1.0 else * frac(i,k) = .0059+.9941*exp(-.003102 * tcel(i,k)*tcel(i,k)) frac(i,k) = .0059+.9941*frac(i,k) endif if (frac(i,k) .lt. 0.01) frac(i,k) = 0. IWP(i,k) = (1.-frac(i,k))*lwcm(i,k)*dp(i,k)*rec_grav*1000. LWP(i,k) = frac(i,k)*lwcm(i,k)*dp(i,k)*rec_grav*1000. end do end do else do k=1,nk do I=1,lmx frac(i,k) = 1. IWP(i,k) = iwcm(i,k)*dp(i,k)*rec_grav*1000. LWP(i,k) = frac(i,k)*lwcm(i,k)*dp(i,k)*rec_grav*1000. end do end do endif * * end do * end do * * * ***************************************************************** * MAIN LOOP FOR MICROPHYSICS SCHEMES ("NEW" OPTICAL PROPERTIES) * ----------------------------------------------------------------- * if (IOPTIX.EQ.2) then * do k=1,nk do I=1,lmx * EQUIVALENT SIZE of PARTICLES * *-------------> determines equivalent size of WATER particles * set number of drops per cm**3 100 for water * and 500 for land if (mg(i) .le. 0.5 .and. ml(i) .le. 0.5) then * cdd(i,k) = 100. rec_cdd(i,k) = 1. / 100. else * cdd(i,k) = 500. rec_cdd(i,k) = 1. / 500. endif * * aird is air density in kg/m3 * aird(i,k) = sig(i,k)*ps(i)/T(i,k)*REC_RGASD * end do end do * CALL VSPOWN1 (REW,(1.+LWCM*1.e4)*LWCM*frac*aird*rec_cdd, + third, nk*lmx) * do k=1,nk do I=1,lmx * REW(i,k) = 3000. * ( (1.+LWCM(i,k)*1.e4)* * + LWCM(i,k)*frac(i,k)*aird(i,k)*rec_cdd(i,k))**third REW(i,k) = 3000. * REW(i,k) REW(i,k) = max( 2.5,REW(i,k)) REW(i,k) = min( 60., REW(i,k)) * * determines array index for given REW ire(i,k) = 2 if ( rew(i,k) .lt. 12.5 ) ire(i,k) = 1 if ( rew(i,k) .gt. 30.0 ) ire(i,k) = 3 * avant d'appeler vspownn il faut definir kwf_ire(i,k) = kwf(2,ire(i,k)) kvf_ire(i,k) = kvf(2,ire(i,k)) ssf_ire(i,k) = ssf(2,ire(i,k)) gwf_ire(i,k) = gwf(2,ire(i,k)) * end do end do * CALL VSPOWNN (kw,rew,kwf_ire,nk*lmx) * do k=1,nk do I=1,lmx * * water diffusivity after Hu and Stammes, JCAM 1993 * * kw = ( kwf(1,ire(i,k))*( rew(i,k)**kwf(2,ire(i,k)) ) + kwf(3,ire(i,k)) )*1.e-3 kw(i,k) = ( kwf(1,ire(i,k)) * kw(i,k) + kwf(3,ire(i,k)) )*1.e-3 * end do end do REI = 15. REC_REI = 1. / 15. KI = .0003 + 1.290 * REC_REI CALL VSEXP (EI,-elsa*ki*IWP,nk*lmx) CALL VSEXP (ENEB,-elsa*ki*IWP-kw*LWP,nk*lmx) * * on elimine une boucle en faisant ceci plus tot CALL VSPOWNN (omw,rew,kvf_ire,nk*lmx) CALL VSPOWNN (ssaw,rew,ssf_ire,nk*lmx) CALL VSPOWNN (gw,rew,gwf_ire,nk*lmx) SSAI = 1.0 - 1.295E-2 - 1.321E-4 *REI * * do k=1,nk do I=1,lmx * * * emissivity of ice after Ebert and Curry, JGR-1992,p3833 * using 11.1 micron parameters as representative * of entire spectrum assume equivalent ice particles * radius of 25 micron will affect both VIS and * IR radiation * * REI = 15. * KI = .0003 + 1.290 / REI * EI = 1. -exp(-elsa*ki *IWP(i,k)) EI(i,k) = 1. - EI(i,k) * * * compute combined ice/water cloud emissivity assuming * the transmission is the product of the ice and water * phase transmissions * * cloud amount is considered outside of the current * main loop due to potential optical depth correction * * cloud emissivity temporarly in ENEB * ENEB(i,k) = 1. - exp( - elsa*ki *IWP(i,k) - kw(i,k) * LWP(i,k) ) ENEB(i,k) = 1. - ENEB(i,k) neb(i,k) = cloud(i,k) * * end do * end do * * * do k=1,nk * do I=1,lmx * OPTICAL THICKNESS * * water optical thickness: Hu and Stammes, JCAM 1993 * for .719 um * * 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 omw(i,k)=LWP(i,k) * ( kvf(1,ire(i,k))*( omw(i,k) ) + kvf(3,ire(i,k)) )*1.e-3 omw(i,k)=omw(i,k)*tuneopw OMI = min(IWP(i,k)*(3.448E-3+2.431*REC_REI) * tuneopi, 25.) OPDEPTH(i,k)= max(omw(i,k) + omi,1.e-10) * * save integrated quantities for output purposes * tlwp (i) = tlwp (i) + lwp(i,k) tiwp (i) = tiwp (i) + iwp(i,k) topthw(i) = topthw(i) + omw(i,k) topthi(i) = topthi(i) + omi * * * SINGLE SCATTERING ALBEDO * * note that this parameter is very close to one for * most of VIS but lowers in near infrared; proper * spectral weighting is important because small changes * will affect substantially the solar heating outgoing * radiance or planetary albedo; It has a smaller influence * on the surface flux. A SSA of unity will create division * by zero in solar code. * SSAI = 1.0 - 1.295E-2 - 1.321E-4 *REI * * WATER Single scattering Albedo: Hu and Stammes, JCAM 1993 * for .719 um * SSAW = 1.- ( ssf(1,ire(i,k))*( rew(i,k)**ssf(2,ire(i,k)) ) + ssf(3,ire(i,k)) ) SSAW(i,k)= 1.- ( ssf(1,ire(i,k))*( SSAW(i,k) ) + ssf(3,ire(i,k)) ) * * weighting by optical depth * OMEGAV(i,k)= (OMW(i,k) * SSAW(i,k) + OMI* SSAI)/ (OMW(i,k)+OMI+eps) OMEGAV(i,k)=max(OMEGAV(i,k),0.9850) OMEGAV(i,k)=min(OMEGAV(i,k),0.999999) * * Ice Asymmetry factor: Ebert and Curry JGR 1992 Eq.8 * for 25 micron particles and a weighting for the * five bands given in their Table 2 * GI = min(0.777 + 5.897E-4 * REI, 0.9999) * * * Water Asymetry factor: Hu and Stammes, JCAM 1993 * for .719 um * GW = ( gwf(1,ire(i,k))*( rew(i,k)**gwf(2,ire(i,k)) ) + gwf(3,ire(i,k)) ) GW(i,k) = ( gwf(1,ire(i,k))*( GW(i,k) ) + gwf(3,ire(i,k)) ) * * weighting by SSA * opt depth * ASYMG(i,k) = (SSAI*OMI*GI + SSAW(i,k)*OMW(i,k)*GW(i,k))/(SSAI*OMI+SSAW(i,k)*OMW(i,k)+eps) * asymg(i,k) = max(asymg(i,k), 0.75 ) * asymg(i,k) = min(asymg(i,k),0.9999) * * geometrical thickness * dz(i,k)= dp(i,k)/aird(i,k)*rec_grav * * end do end do else if (IOPTIX.EQ.1) then * * **************************************************************** * MAIN LOOP FOR CONVENTIONAL CONDENSATION SCHEMES * ("OLD" OPTICAL PROPERTIES) * ---------------------------------------------------------------- * REI = 15. REC_REI = 1. / 15. KI = .0003 + 1.290 * REC_REI * CALL VSEXP (EW,-0.087*elsa*LWP,nk*lmx) CALL VSEXP (EI,-elsa*ki*IWP,nk*lmx) * do k=1,nk do I=1,lmx * * emissivity of water after Stephens, JAS 78 * we take a 0.144 factor as average of his * downward and upward emissivity which divided * by diffusivity factor elsa yields 0.087 * * EW = 1. -exp(-0.087*elsa* LWP(i,k)) EW(i,k) = 1. - EW(i,k) * * emissivity of ice after Ebert and Curry, JGR-1992,p3833 * using 11.1 micron parameters as representative * of entire spectrum assume equivalent ice particles * radius of 25 micron will affect both VIS and * IR radiation * * EI(i,k) = 1. -exp(-elsa*ki *IWP(i,k)) EI(i,k) = 1. - EI(i,k) * * compute combined ice/water cloud emissivity assuming * the transmission is the product of the ice and water * phase transmissions * EMISS = 1. - (1.-EI(i,k))* (1.-EW(i,k)) * * effective cloud cover * ENEB(i,k)= cloud(i,k)*emiss * * black clouds for istcond non 3 * * eneb(i,k)= cvmgt(cloud(i,k),eneb(i,k),istcond.NE.3) * * optical thickness at nadir is computed * set number of drops per cm**3 100 for water * and 500 for land * * The following line is commented out to ensure bitwise * validation of the GEMDM global model. It should be * activated in a future version of the code. if (mg(i).le.0.5) then * cdd(i,k) = 100. rec_cdd(i,k) = 1. / 100. else * cdd(i,k) = 500. rec_cdd(i,k) = 1. / 500. endif * * aird is air density in kg/m3 aird(i,k) = sig(i,k)*ps(i)/T(i,k)*REC_RGASD * end do end do * CALL VSPOWN1 (REW,LWCM*frac*aird*rec_cdd,third,nk*lmx) * do k=1,nk do I=1,lmx * * this parameterization from H. Barker, based * on aircraft data range 4-17 micron is that specified * by Slingo for parameterizations * * REW(i,k) = max(4., 754.6 * (LWCM(i,k)*frac(i,k)*aird(i,k)*rec_cdd(i,k))**third) REW(i,k) = max(4., 754.6 * REW(i,k)) REW(i,k) = min(REW(i,k),17.) * * slingo JAS 89 p 1420 for weighted average of * bands 1-5 (all spectrum) * OMW(i,k) = LWP(i,k)* (2.622E-2 + 1.356/REW(i,k) ) * * follows Ebert and Curry, 1992 * no variation as function of band * ice optical depth limited to 25; water to 125. * omw(i,k)= min(omw(i,k)*tuneopw,125.) OMI = min(IWP(i,k)*(3.448E-3+2.431*REC_REI) * tuneopi, 25.) OPDEPTH(i,k)= max(omw(i,k) + omi,1.e-10) * * * save integrated quantities for output purposes * tlwp (i) = tlwp (i) + lwp(i,k) tiwp (i) = tiwp (i) + iwp(i,k) topthw(i) = topthw(i) + omw(i,k) topthi(i) = topthi(i) + omi * SSAI = 1.0 - 1.295E-2 - 1.321E-4 *REI * Slingo JAS 1989 for weighted average bands 1-4 p 1420 SSAW(i,k) = 1.0 - 6.814E-3 - 4.842E-4 *REW(i,k) * * * weighting by optical depth * OMEGAV(i,k)= (OMW(i,k) * SSAW(i,k) + OMI* SSAI)/ (OMW(i,k)+OMI+eps) OMEGAV(i,k)=max(OMEGAV(i,k),0.9850) OMEGAV(i,k)=min(OMEGAV(i,k),0.9999) * * * Ice Asymmetry factor: Ebert and Curry JGR 1992 Eq.8 * for 25 micron particles and a weighting for the * five bands given in their Table 2 * GI = min(0.777 + 5.897E-4 * REI, 0.9999) * * Water Asymetry factor: Slingo 1989 * GW(i,k) = min(0.804 + 3.850E-3*REW(i,k), 0.9999) * * weighting by SSA * opt depth * ASYMG(i,k) = (SSAI*OMI*GI + SSAW(i,k)*OMW(i,k)*GW(i,k))/(SSAI*OMI+SSAW(i,k)*OMW(i,k)+eps) * asymg(i,k) = max(asymg(i,k),GI) * asymg(i,k) = min(asymg(i,k),0.9999) * * geometrical thickness * dz(i,k)= dp(i,k)/aird(i,k)*rec_grav * end do end do * * endif * * * * ************************************************************ * END OF MAIN LOOPS * ----------------------------------------------------------- * if (IOPTIX.EQ.2) then * * ****************************************************************** * CORRECTION FOR MAXIMUM TOTAL OPTICAL DEPTH & * FINAL CLOUD*EMISSIVITY CALCULATIONS * ------------------------------------------------------------------ * * temporary use of ipbl as a work field. * ipbl is 1. when there is no correction * otherwise ipbl is a number between zero and one * do i=1, lmx ipbl(i)=min( 1., 20./( max(topthi(i)+topthw(i), 1.e-10) ) ) enddo * only optical depth and cloud emissivity * need to be re-scaled * do k=1, nk do i=1, lmx OPDEPTH(i,k) = ipbl(i) * OPDEPTH(i,k) eneb(i,k) = neb(i,k) * $ ( 1. - ( 1.-eneb(i,k) )**ipbl(i) ) enddo enddo * * endif * * Diagnostics : cloud top pressure (ctp) and temperature (ctt) * do i=1,lmx ctp (i) = 110000. ctt (i) = 310. top(i) = .true. transmissint(i,1) = 1. - eneb(i,1) if ( (1.-transmissint(i,1)) .gt. 0.99 .and. top(i) ) then ctp(i) = sig(i,1)*ps(i) ctt(i) = t(i,1) top(i) = .false. end if end do do k=2,nk do i=1,lmx transmissint(i,k) = transmissint(i,k-1) * (1.-eneb(i,k)) if ( (1.-transmissint(i,k)) .gt. 0.99 .and. top(i) ) then ctp(i) = sig(i,k)*ps(i) ctt(i) = t(i,k) top(i) = .false. end if end do end do * * * ****************************************************************** * AEROSOLS LOADING * ------------------------------------------------------------------ * do 10 I =1,lmx sdz(i) = 0. ipbl(i) = 0. 10 continue * do 5 k=nk,1,-1 do 6 I=1,lmx if ( int(ipbl(i)).eq.0 ) then * pbl heiht recomputed as sum of layer thicknesses sdz(i)=sdz(i)+dz(i,k) * level closest to pbl if (sdz(i) .gt. pbl(i)) ipbl(i) = float(k) endif 6 continue 5 continue * * distributing aerosols * optical thickness higher over land; * decreases with latitude * See Toon and Pollack, J. Appl. Meteor, 1976, p.235 * aero being total optical thickness, we distribute it * within PBL in proportion with geometrical thickness * above pbl aerosol optical depth is kept negligible * ct=2./pi REC_180=1./180. do 3 k=1,nk do 4 i=1,lmx tauae(i,k,1)=1.e-10 tauae(i,k,2)=1.e-10 rlat = lat(i)*pi*REC_180 if ( k.ge.INT(ipbl(i)) ) then * if ( mg(i).ge.0.5.or.ml(i).ge.0.5) then * * over land * aero = 0.25 - 0.2*ct * abs(lat(i)) tauae(i,k,1)= max(aero * dz(i,k)/sdz(i), 1.e-10) else * * over ocean * aero = 0.13 - 0.1*ct * abs(lat(i)) tauae(i,k,2)= max(aero *dz(i,k)/sdz(i), 1.e-10) endif endif * * other types of aerosols set to negligible * tauae(i,k,3)=1.e-10 tauae(i,k,4)=1.e-10 tauae(i,k,5)=1.e-10 4 continue 3 continue return end @ 4.1 log @bidon @ text @d60 2 d863 1 a863 10 if ( (ioptix.eq.2 .and. + (mg(i).ge.0.5.or.ml(i).ge.0.5)) .or. + (ioptix.lt.2 .and. + mg(i).ge.0.5 ) ) then * * The following line is commented out to ensure bitwise * validation of the GEMDM global model. It should be * activated in a future version of the code, and * the previous lines should be deleted. c if ( mg(i).ge.0.5.or.ml(i).ge.0.5) then @ 4.0 log @La version 4.0 de la physique a ete creee le 4 septembre 2003. C'est la premiere version de la programmatheque des parametrages destinee a etre mise en exploitation sur le nouveau calculateur IBM. Principaux changements : ---------------------- * Correction de bogues introduits lors de la conversion vers IBM * Nouvelles optimisations pour IBM * Modifications pour permettre l'interpolation des profils atmospheriques pres de la surface @ text @@ 3.9 log @La version 3.9 de la physique a ete creee le 16 juin 2003. Elle constitue la premiere version de conversion vers le calculateur IBM. Le nouveau code de "gravity wave drag" sgoflx3.ftn est une copie du code linearise lin_sgoflx1.ftn. @ text @d59 1 d458 2 a459 2 CALL EXPONEN4 (REW,(1.+LWCM*1.e4)*LWCM*frac*aird*rec_cdd, + third, nk*lmx, nk*lmx, 1) d474 1 a474 1 * pour faire fonctionner exponen4 il faut definir d483 1 a483 1 CALL EXPONEN4 (kw,rew,kwf_ire,nk*lmx,nk*lmx,nk*lmx) d503 3 a505 3 CALL EXPONEN4 (omw,rew,kvf_ire,nk*lmx,nk*lmx,nk*lmx) CALL EXPONEN4 (ssaw,rew,ssf_ire,nk*lmx,nk*lmx,nk*lmx) CALL EXPONEN4 (gw,rew,gwf_ire,nk*lmx,nk*lmx,nk*lmx) a540 4 * CALL EXPONEN4 (omw,rew,kvf_ire,nk*lmx,nk*lmx,nk*lmx) * CALL EXPONEN4 (ssaw,rew,ssf_ire,nk*lmx,nk*lmx,nk*lmx) * CALL EXPONEN4 (gw,rew,gwf_ire,nk*lmx,nk*lmx,nk*lmx) * SSAI = 1.0 - 1.295E-2 - 1.321E-4 *REI d688 1 a688 1 CALL EXPONEN4 (REW,LWCM*frac*aird*rec_cdd,third,nk*lmx,nk*lmx,1) a786 2 CALL EXPONEN4 (eneb,1.-eneb,ipbl,nk*lmx,nk*lmx,lmx) * d790 2 a791 2 * eneb(i,k) = neb(i,k) * ( 1. - ( 1.-eneb(i,k) )**ipbl(i) ) eneb(i,k) = neb(i,k) * ( 1. - eneb(i,k) ) @ 3.8 log @La version 3.8 de la physique a ete creee le 12 mars 2003 Principaux changements par rapport a la version 3.72 : ---------------------------------------------------- * contenu de la rustine r2 de la version 3.72 * code developpe pour le modele regional a 15 km - MOISTKE (refonte) - MIXPHASE (avec BLG) - KTRSNT - proprietes optiques des nuages * option ADVECTKE reactivee * BOUJO disponible dans eturbl * ajouts importants au code de physique linearisee * nouvelles cles : AS,BETA,KKL,KFCPCP,SHLCVT(2),STRATOS,QCO2 * nombreux diagnostics supplementaires * optimisation des series temporelles * diagnostics supplementaires pour CHRONOS et AURAMS * correction d'une multitude de bogues mineurs @ text @d54 7 a124 3 #if defined (CVMG) #include "cvmg.cdk" #endif d142 5 a146 4 integer i,k, ire logical top real rew, kw, lwcm1, tuneopw, omw, gw, ssaw real rei, ki, iwcm1, tuneopi, omi, gi, ssai d148 2 a149 2 real ei, ew, elsa, emiss real aird, third, tcel, cdd d152 4 d168 3 a273 1 lwc(i,k)=0.4*asymg(i,k) d277 5 a281 1 lwc(i,k)=cvmgt(0.,lwc(i,k),sig(i,k).lt.0.050) d291 3 a293 1 lwc(i,k)=cvmgt(0.,lwc(i,k),sig(i,k).lt.0.070) d329 5 a333 2 cloud(i,k) = cvmgt(max(cloud(i,k) ,0.01),0.0, + (lwcm(i,k)+iwcm(i,k)).gt.1.e-6) d336 4 a339 2 lwcm(i,k)= cvmgt(0., lwcm(i,k), cloud(i,k) .lt. 0.01) iwcm(i,k)= cvmgt(0., iwcm(i,k), cloud(i,k) .lt. 0.01) d369 8 a376 2 dp(i,k)=cvmgt(dp2,dp1,k.eq.1) dp(i,k)=cvmgt(dp3,dp(i,k),k.eq.nk) d378 5 d397 11 a407 4 tcel=T(i,k)-TCDK frac(i,k) = cvmgt(1.,.0059+.9941*exp(-.003102 * tcel*tcel), x tcel.ge.0.) frac(i,k)= cvmgt(0.,frac(i,k), frac(i,k).lt.0.01) d409 4 a412 1 IWP(i,k) = (1.-frac(i,k))*lwcm(i,k)*dp(i,k)/grav*1000. d414 2 a415 1 * d417 4 a420 1 IWP(i,k) = iwcm(i,k)*dp(i,k)/grav*1000. a422 1 LWP(i,k) = frac(i,k)*lwcm(i,k)*dp(i,k)/grav*1000. d424 2 a425 2 end do end do d442 7 a448 1 cdd =cvmgt(100.,500.,mg(i).le.0.5.and.ml(i).le.0.5) d452 7 a458 1 aird = sig(i,k)*ps(i)/T(i,k)/RGASD d460 7 a466 4 REW = 3000. * ( (1.+LWCM(i,k)*1.e4)* + LWCM(i,k)*frac(i,k)*aird/cdd )**third REW = max( 2.5,REW) REW = min( 60., REW) d470 13 a482 4 ire = 2 if ( rew .lt. 12.5 ) ire = 1 if ( rew .gt. 30.0 ) ire = 3 d484 2 d488 6 d495 16 a510 2 kw = ( kwf(1,ire)*( rew**kwf(2,ire) ) + kwf(3,ire) )*1.e-3 d518 5 a522 3 REI = 15. KI = .0003 + 1.290 / REI EI = 1. -exp(-elsa*ki *IWP(i,k)) a523 1 d533 2 a534 1 ENEB(i,k) = 1. - exp( - elsa*ki *IWP(i,k) - kw * LWP(i,k) ) d537 10 d551 4 a554 3 omw=LWP(i,k) * ( kvf(1,ire)*( rew**kvf(2,ire) ) + kvf(3,ire) )*1.e-3 omw=omw*tuneopw d556 1 a556 1 OMI = min(IWP(i,k)*(3.448E-3+2.431/REI) * tuneopi, 25.) d558 1 a558 1 OPDEPTH(i,k)= max(omw + omi,1.e-10) d564 1 a564 1 topthw(i) = topthw(i) + omw d578 1 a578 1 SSAI = 1.0 - 1.295E-2 - 1.321E-4 *REI d584 2 a585 1 SSAW= 1.- ( ssf(1,ire)*( rew**ssf(2,ire) ) + ssf(3,ire) ) d589 1 a589 1 OMEGAV(i,k)= (OMW * SSAW + OMI* SSAI)/ (OMW+OMI+eps) d603 2 a604 1 GW = ( gwf(1,ire)*( rew**gwf(2,ire) ) + gwf(3,ire) ) d608 1 a608 1 ASYMG(i,k) = (SSAI*OMI*GI + SSAW*OMW*GW)/(SSAI*OMI+SSAW*OMW+eps) d616 1 a616 1 dz(i,k)= dp(i,k)/(aird*grav) d629 7 d644 2 a645 1 EW = 1. -exp(-0.087*elsa* LWP(i,k)) d653 2 a654 3 REI = 15. KI = .0003 + 1.290 / REI EI = 1. -exp(-elsa*ki *IWP(i,k)) d660 1 a660 1 EMISS = 1. - (1.-EI)* (1.-EW) d677 7 a683 2 c cdd =cvmgt(100.,500.,mg(i).le.0.5.and.ml(i).le.0.5) cdd =cvmgt(100.,500.,mg(i).le.0.5) d686 10 a695 2 aird = sig(i,k)*ps(i)/T(i,k)/RGASD d700 3 a702 2 REW = max(4., 754.6 * (LWCM(i,k)*frac(i,k)*aird/cdd)**third) REW = min(REW,17.) d708 1 a708 1 OMW = LWP(i,k)* (2.622E-2 + 1.356/REW ) d714 3 a716 3 omw= min(omw*tuneopw,125.) OMI = min(IWP(i,k)*(3.448E-3+2.431/REI) * tuneopi, 25.) OPDEPTH(i,k)= max(omw + omi,1.e-10) d723 1 a723 1 topthw(i) = topthw(i) + omw d728 1 a728 1 SSAW = 1.0 - 6.814E-3 - 4.842E-4 *REW d733 1 a733 1 OMEGAV(i,k)= (OMW * SSAW + OMI* SSAI)/ (OMW+OMI+eps) d746 1 a746 1 GW = min(0.804 + 3.850E-3*REW, 0.9999) d750 1 a750 1 ASYMG(i,k) = (SSAI*OMI*GI + SSAW*OMW*GW)/(SSAI*OMI+SSAW*OMW+eps) d758 1 a758 1 dz(i,k)= dp(i,k)/(aird*grav) d790 2 d795 2 a796 2 eneb(i,k) = neb(i,k) * $ ( 1. - ( 1.-eneb(i,k) )**ipbl(i) ) d844 1 a844 1 ipbl(i)=cvmgt(float(k),ipbl(i),sdz(i).gt.pbl(i)) d858 1 d863 1 a863 1 rlat = lat(i)*pi/180. @