Warning, /verification/atm_gray/input/calc_Qsat.m is written in an unsupported language. File is not indexed.
view on githubraw file Latest commit 42875574 on 2024-11-26 20:20:16 UTC
4287557481 Jean*0001 function [qsat,dQdT]=calc_Qsat(kpot,p,t);
0002 % [qsat,dQdT]=calc_Qsat(kpot,p,t);
0003 % kpot=0 : t = T = absolute Temp ; kpot=1 : t = theta = pot.temp
0004 % p = normalized pressure = p/1000.mb
0005 % output: qsat (g/kg) ; d.qsat/d.T
0006 % ktyp=1 : p in mb ; t = absolute Temp ;
0007
0008 g=9.81;
0009 kappa=2./7.;
0010
0011 np=prod(size(p));
0012 nt=prod(size(t));
0013 nx=nt/np;
0014 if nx ~= fix(nx) | nx < 1
0015 fprintf(' Pb in dimensions: dim_p,dim_t= %i , %i ; ratio= %f \n',np,nt,nx);
0016 qsat=0.; dQdT=0.;
0017 return
0018 end
0019
0020 T=reshape(t,nx,np);
0021 P=reshape(p,1,np);
0022 if kpot == 1,
0023 factp=P.^kappa;
0024 fp=ones(nx,1)*factp;
0025 size(fp);
0026 T=fp.*T;
0027 end
0028
0029 e0=6.108e-3 ; c1=17.269 ; c2=21.875 ; t0c=273.16 ; ct1=35.86 ; ct2=7.66 ;
0030 c0q=622. ; c1q=0.378 ;
0031
0032 qsat=zeros(nx,np);
0033 dQdT=zeros(nx,np);
0034
0035 for k=1:np,
0036 for j=1:nx,
0037 if T(j,k) >= t0c
0038 qE=e0*exp(c1*(T(j,k)-t0c)/(T(j,k)-ct1)) ;
0039 qsat(j,k)=c0q*qE/(P(k)-c1q*qE);
0040 dQdT(j,k)=c1*P(k)*qsat(j,k)/(P(k)-c1q*qE)*(t0c-ct1)/(T(j,k)-ct1)^2 ;
0041 elseif T(j,k) > ct2
0042 qE=e0*exp(c2*(T(j,k)-t0c)/(T(j,k)-ct2)) ;
0043 qsat(j,k)=c0q*qE/(P(k)-c1q*qE);
0044 dQdT(j,k)=c2*P(k)*qsat(j,k)/(P(k)-c1q*qE)*(t0c-ct2)/(T(j,k)-ct2)^2 ;
0045 end
0046 end
0047 end
0048
0049 return
