program buddiag c... This module buddiag calculates the large scale budgets and c related diagnostics from cloud model outputs. The I/O and c calculations are for the 2d Szeto-Cho model results and c therefore substantial modification to the code will be needed c in order to apply it to the outputs of other models. This c version assumes that the results have reached a quasi-steady c state and that temporal averaging can be neglected. c c revision: Aug, 27, 95. parameter(m1=513, m2=1, m3=41) common/bk1/ qc(m1,m3), qi(m1,m3), qr(m1,m3), qs(m1,m3), & qg(m1,m3), qv(m1,m3) common/bk2/ t(m1,m3), u(m1,m3), v(m1,m3), w(m1,m3), & th(m1,m3) common/bk3/ dx,dy,dz,dt,dxx,dyy,dzz,dtt,rho(m3), pe(m3), & pcpn(m1) common/bk4/ is,ie,js,je,ks,ke common/bk5/ tm(m3),um(m3),vm(m3),wm(m3),qim(m3),qcm(m3), & qrm(m3),qsm(m3),qgm(m3),qvm(m3),iwcm(m3),lwcm(m3) common/bk6/ td(m1,m3), ud(m1,m3), vd(m1,m3), wd(m1,m3), & thd(m1,m3) common/bk7/ tv(m3),uv(m3),vv(m3),wv(m3),qiv(m3),qcv(m3), & qrv(m3),qsv(m3),qgv(m3),qvv(m3),iwcv(m3),lwcv(m3) common/bk8/ lwpm(m3),iwpm(m3),lwpv(m3),iwpv(m3),icloud(m1) real lwcm,lwpm,lwcv,lwpv,iwcm,iwpm,iwcv,iwpv real z(m3),zz(m3),q1(m3),q2(m3),qu(m3),qqv(m3),rhom(m3) real sigma(m3),sigmac(m3),taucm(m3),taucv(m3) real tauim(m3),tauiv(m3),thm(m3),rvm(m3),qdot(m1,m3) real tauc(m1,m3),taui(m1,m3),taul(m1),qdot2(m1,m3) real ueddy(m3),veddy(m3),weddy(m3),theddy(m3),rveddy(m3) real sw1(m3),sw2(m3) ********** VARIABLES *************** c...m1,m2,m3 are the total number of grids in the 3 directions c of the CRM model domain c qc, qi, qr, qs, qg, qv are mixing ratios of cloud water, c cloud ice, rain, snow, graupel and water vapor respectively. c lwc,iwc= liquid water content, ice water content c t = temperature c th = potential temperature c u,v,w = the 3 velocity components in a Catesian co-ordinate c qdot, qdot2 = total latent heating and total latent heating c w/o fusion c pcpn = surface precip rate c rho = mean air density and c pi = mean non-dimensional pressure (defined as Poo/P)**(R/Cp) c here) c dx,dy,dz,dt = grid intervals and timestep in CRM c dxx, dyy, dzz,dtt defines the subdomain szie and time period c for the diagnostic calculations (see DX,DY,DZ,DT in the following) c Xm = mean of varibale X c Xd = fluctuations of X c Xv = variance of X *************************************************************** c c LARGE SCALE IMPACT CALCULATIONS: c c Revision: July, 14, 1995 c c A. Cloud-Resolving Model (CRM) outputs that are useful to c GCM modellers: c c Over subregions of the CRM with dimensions similar to those c of a typical GCM grid box (~300x300 sq. km x 500 m) and over c a time interval representative of the typical time step used c in a GCM integration, the following statistics (as function c of height)can be generated from the CRM outputs: c c (1) Mean and variance of temperature c c (2) Mean and variance of mixing ratios of various water c substance, layer liquid water content, layer liquid water c path, layer optical depth and surface precipitation rate c c (3) Cloud coverage statistics c c (4) Mean long (short) wave fluxes c c (5) Fluctuations in the velocity fields and associated eddy c fluxes of heat, moisture and momentum c c (6) Apparent heat, moisture and momentum sources/sinks c c Calculations over subdomains with and without embedded c convection (and/or CSI) would be useful. It is also c understood that (2) and (3) will be calculated for the c model-resolved cloud types only. c c B. Calculations of the deliverables: c c Let dx', dy', dz' and dt' be the grid-lengths and time step c used in the CRM and DX, DY, DZ, DT be the typical grid- c lengths and time step used in a GCM (DX, DY ~ 300 km, DZ ~ c 500m, DT ~ 30 min), then the mean value of the CRM c prognostic variable p can be defined by: c c mean (p) = sum(p dx'dy'dz'dt') / (DXDYDZDT) c c where "sum" represents summation over all the grid points c within the subdomain DX.DY.DZ. of the CRM grid box and over c a time interval DT. Likewise, the variance of p can be c defined by: c c var (p) = sum( (p-mean (p))**2 dx'dy'dz'dt') / (DXDYDZDT) c c c . The optical depth tau is defined as: c c tau = INTh INTr [f(r) Qe(s,n(l) pi r**2 ] c c where "INTx" represents an integral over x. The first c integral is over the cloud depth h and the second integral c is over the cloud droplet radius r. f(r) is the spectral c density of droplets of radius r while Qe is the extinction c factor (ratio of extinction to the cross-sectional area of c cloud droplets) as a function of the size parameter s and c refractive index n at wavelength l. c c The functional form of f and Qe typically varies from CRM to c CRM. c c . The liquid water path W is defined as the integrated c liquid water through a cloud of depth h: c c W = INTh [rho rc] c c where rho is air density and rc is the cloud water mixing c ratio. c c . The cloud coverage fraction can be calculated if we c specify the threshold value of cloud water (plus cloud ice) c mixing ratio to be 0.001 g kg^-1. This should allow for c cirrus. The cloud fraction (sigma) can then be calculated c by: c c sigma = sum(dx'dy') / (DX DY) c c where the summation is over all grid points where the cloud c mixing ratios exceeds the respective threshold mixing c ratios. c c . It would be of interest to calculate the fractional c coverage of convective clouds if moist convection is c resolved explicitly in the model. This can be done if we c specify a threshold convective cloud water (plus cloud ice) c mixing ratio. There is no precise value for this threshold c and it will undoubtedly vary with case anyway. For many of c our Canadian cases, a value of say 0.2 g kg^-1 is reasonable c for instance. Another condition that might be added is that c the in-cloud updraft speed be greater than 0.4 m/s. c c . For the purpose of comparing model cloud coverage to c satellite data (eg. ISCCP), the cloud coverage can also be defined c as the area where the column visible optical thickness be c greater than 0.1. c c . The fluctuations for velocity component i are defined by: c c ui" = ui - mean (ui) c c . The associated eddy fluxes of momentum, heat and moisture c are: c c mean (ui" uj") c mean (ui" T") c mean (ui" q") c c where <> denotes horizontal averages over the model domain c and T is potential temperature while q is the water vapor c mixing ratio. c c . The apparent heat and moisture sources and sinks are c defined as: c c Q1 = mean (DT/Dt) = - / mean(@ (w"T")) + Q c + DQ1 c c Q2 = -/Cp mean (Dq/Dt) = /Cp mean(@ (w"q")) c + Q + DQ2 c c where E is Exner pressure ( Cp (P/Poo)**(Rd/Cp)), Cp = heat c capacity of dry air, P = pressure and Poo = reference c pressure, Rd = gas constant for dry air. L = latent heat of c vaporization, @() means partial derivative with respective c to the vertical co-ordinate z, Q = latent heating and DQi c denotes subgrid scale diffusions. c c . The apparent sources and sinks of horizontal momenta can c be defined as: c c Qui = mean(Dui/Dt) = - (1/) mean(@( w"ui")) + Dqui c ********************************************************************* c.. Inputs: print*,'Define the subdomain for the diagnostic calculations:' print*,'i,j,k are grid indexes in the 3 directions:' print*,'Enter the starting and ending i grid indexes: ' read*,is,ie print*,'Enter the starting and ending j grid indexes: ' read*,js,je print*,'Enter the starting and ending k grid indexes: ' read*,ks,ke c..Assign constants .. dx = 15000.0 dy = 0.0 dz = 250.0 dy = 250.0 c dxx = 300000.0 dxx = float(ie-is)*dx dyy = 0.0 dzz = 500.0 dtt = 1200.0 kl = int(dzz/dz) tot = float((ie-is+1)*(je-js+1)*(kl+1)) do k=ks,ke z(k)=float(k-1)*dz/1000. enddo kk=0 do k=ks,ke,kl kk=kk+1 zz(kk)=(float(k+k+kl)*0.5-1.)*dz/1000. enddo c... Read in CRM outputs ... open(2, form='unformatted', file='THTA01') do k=1,m3 do i=1,m1 read(2) th(i,k) enddo enddo close(2) open(30,file='THTATT',form='unformatted') do k=1,m3 read(30) dumt do i=1,m1 th(i,j)=th(i,j)+dumt enddo enddo close(30) open(2, form='unformatted', file='TEMP01') do k=1,m3 do i=1,m1 read(2) t(i,k) enddo enddo close(2) open(2, form='unformatted', file='XVEL01') do k=1,m3 do i=1,m1 read(2) u(i,k) enddo enddo close(2) open(2, form='unformatted', file='YVEL01') do k=1,m3 do i=1,m1 read(2) v(i,k) enddo enddo close(2) open(2, form='unformatted', file='ZVEL01') do k=1,m3 do i=1,m1 read(2) w(i,k) enddo enddo close(2) open(2, form='unformatted', file='VMIX01') do k=1,m3 do i=1,m1 read(2) qv(i,k) enddo enddo close(2) open(2, form='unformatted', file='CMIX01') do k=1,m3 do i=1,m1 read(2) qc(i,k) enddo enddo close(2) open(2, form='unformatted', file='ZMIX01') do k=1,m3 do i=1,m1 read(2) qi(i,k) enddo enddo close(2) open(2, form='unformatted', file='SMIX01') do k=1,m3 do i=1,m1 read(2) qs(i,k) enddo enddo close(2) open(2, form='unformatted', file='GMIX01') do k=1,m3 do i=1,m1 read(2) qg(i,k) enddo enddo close(2) open(2, form='unformatted', file='RMIX01') do k=1,m3 do i=1,m1 read(2) qr(i,k) enddo enddo close(2) open(2, form='unformatted', file='QDOT01') do k=1,m3 do i=1,m1 read(2) qdot(i,k) enddo enddo close(2) open(2, form='unformatted', file='QDOT02') do k=1,m3 do i=1,m1 read(2) qdot2(i,k) enddo enddo close(2) open(2, file='SFRF01') do i=1,m1 read(2,101) pcpn(i) enddo close(2) 100 format(e16.9) 101 format(e16.9,e16.9) open(33,file='INICO') do k=1,m3 read(33,100) dum enddo do k=1,m3 read(33,100) dum enddo do k=1,m3 read(33,100) rho(k) enddo do k=1,m3 read(33,100) pe(k) enddo close(33) C.... Calculate mean and variance of temperature: kkt=0 kk=0 do k=ks,ke,kl kk=kk+1 tm(kk)=0.0 do l=k,k+kl do i=is,ie tm(kk)=t(i,l)+tm(kk) enddo enddo tm(kk)=tm(kk)/tot kkt=kkt+1 enddo kk=0 do k=ks,ke,kl kk=kk+1 tv(kk)=0.0 do l=k,k+kl do i=is,ie tv(kk)=(t(i,l)-tm(kk))**2 + tv(kk) enddo enddo tv(kk)=tv(kk)/tot enddo open(3,file='tempm') do k=1,kkt write(3,101) tm(k),zz(k) enddo close(3) open(3,file='tempv') do k=1,kkt write(3,101) tv(k),zz(k) enddo close(3) C.... Calculate mean and variance of cloud water mixing ratio kk=0 do k=ks,ke,kl kk=kk+1 qcm(kk)=0.0 rhom(kk)=0. do l=k,k+kl rhom(kk)=rho(k)+rhom(kk) do i=is,ie qcm(kk)=qc(i,l)+qcm(kk) enddo enddo rhom(kk)=rhom(kk)/float(kl+1) qcm(kk)=qcm(kk)/tot lwcm(kk)=rhom(kk)*qcm(kk) enddo kk=0 do k=ks,ke,kl kk=kk+1 qcv(kk)=0.0 do l=k,k+kl do i=is,ie qcv(kk)=(qc(i,l)-qcm(kk))**2 + qcv(kk) enddo enddo qcv(kk)=qcv(kk)/tot lwcv(kk)=rhom(kk)*qcv(kk) enddo open(3,file='cloudm') do k=1,kkt write(3,101) qcm(k),zz(k) enddo close(3) open(3,file='cloudv') do k=1,kkt write(3,101) qcv(k),zz(k) enddo close(3) open(3,file='lwcm') do k=1,kkt write(3,101) lwcm(k),zz(k) enddo close(3) open(3,file='lwcv') do k=1,kkt write(3,101) lwcv(k),zz(k) enddo close(3) C.... Calculate mean and variance of cloud ice mixing ratio kk=0 do k=ks,ke,kl kk=kk+1 qim(kk)=0.0 do l=k,k+kl do i=is,ie qim(kk)=qi(i,l)+qim(kk) enddo enddo qim(kk)=qim(kk)/tot iwcm(kk)=rhom(kk)*qim(kk) enddo kk=0 do k=ks,ke,kl kk=kk+1 qiv(kk)=0.0 do l=k,k+kl do i=is,ie qiv(kk)=(qi(i,l)-qim(kk))**2 + qiv(kk) enddo enddo qiv(kk)=qiv(kk)/tot iwcv(kk)=rhom(kk)*qiv(kk) enddo open(3,file='icem') do k=1,kkt write(3,101) qim(k),zz(k) enddo close(3) open(3,file='icev') do k=1,kkt write(3,101) qiv(k),zz(k) enddo close(3) open(3,file='iwcm') do k=1,kkt write(3,101) iwcm(k),zz(k) enddo close(3) open(3,file='iwcv') do k=1,kkt write(3,101) iwcv(k),zz(k) enddo close(3) C.... Calculate mean and variance of rain mixing ratio kk=0 do k=ks,ke,kl kk=kk+1 qrm(kk)=0.0 do l=k,k+kl do i=is,ie qrm(kk)=qr(i,l)+qrm(kk) enddo enddo qrm(kk)=qrm(kk)/tot enddo kk=0 do k=ks,ke,kl kk=kk+1 qrv(kk)=0.0 do l=k,k+kl do i=is,ie qrv(kk)=(qr(i,l)-qrm(kk))**2 + qrv(kk) enddo enddo qrv(kk)=qrv(kk)/tot enddo open(3,file='rainm') do k=1,kkt write(3,101) qrm(k),zz(k) enddo close(3) open(3,file='rainv') do k=1,kkt write(3,101) qrv(k),zz(k) enddo close(3) C.... Calculate mean and variance of snow mixing ratio kk=0 do k=ks,ke,kl kk=kk+1 qsm(kk)=0.0 do l=k,k+kl do i=is,ie qsm(kk)=qs(i,l)+qsm(kk) enddo qsm(kk)=qsm(kk)/tot enddo enddo kk=0 do k=ks,ke,kl kk=kk+1 qsv(kk)=0.0 do l=k,k+kl do i=is,ie qsv(kk)=(qs(i,l)-qsm(kk))**2 + qsv(kk) enddo qsv(kk)=qsv(kk)/tot enddo enddo kk=0 open(3,file='snowm') do k=1,kkt write(3,101) qsm(k),zz(k) enddo close(3) open(3,file='snowv') do k=1,kkt write(3,101) qsv(k),zz(k) enddo close(3) C.... Calculate mean and variance of graupel mixing ratio do k=ks,ke,kl kk=kk+1 qgm(kk)=0.0 do l=k,k+kl do i=is,ie qgm(kk)=qg(i,l)+qgm(kk) enddo enddo qgm(kk)=qgm(kk)/tot enddo kk=0 do k=ks,ke,kl kk=kk+1 qgv(kk)=0.0 do l=k,k+kl do i=is,ie qgv(kk)=(qg(i,l)-qgm(kk))**2 + qgv(kk) enddo enddo qgv(kk)=qgv(kk)/tot enddo open(3,file='graupelm') do k=1,kkt write(3,101) qgm(k),zz(k) enddo close(3) open(3,file='graupelv') do k=1,kkt write(3,101) qgv(k),zz(k) enddo close(3) C.... Calculate mean and variance of water vapor mixing ratio kk=0 do k=ks,ke,kl kk=kk+1 qvm(kk)=0.0 do l=k,k+kl do i=is,ie qvm(kk)=qv(i,l)+qvm(kk) enddo enddo qvm(kk)=qvm(kk)/tot enddo kk=0 do k=ks,ke,kl kk=kk+1 qvv(kk)=0.0 do l=k,k+kl do i=is,ie qvv(kk)=(qv(i,l)-qvm(kk))**2 + qvv(kk) enddo enddo qvv(kk)=qvv(kk)/tot enddo open(3,file='vaporm') do k=1,kkt write(3,101) qvm(k),zz(k) enddo close(3) open(3,file='vaporv') do k=1,kkt write(3,101) qvv(k),zz(k) enddo close(3) C.... Calculate mean and variance of surface precip rate pcpnm=0.0 do i=is,ie pcpnm=pcpn(i)+pcpnm enddo pcpnm=pcpnm/float(ie-is+1) pcpnv=0.0 do i=is,ie pcpnv=(pcpn(i)-pcpnm)**2 + pcpnv enddo pcpnv=pcpnv/float(ie-is+1) open(2,file='diag.out') write(2,201) pcpnm 201 format('mean surface precip rate = ',e16.9) write(2,202) pcpnv 202 format('variance of surface precip rate = ',e16.9) c.... Calculate layer liquid (ice) water path statistics do k=1,kkt lwpm(k) = lwcm(k) * float(kl)*dz lwpv(k) = lwcv(k) * float(kl)*dz iwpm(k) = iwcm(k) * float(kl)*dz iwpv(k) = iwcv(k) * float(kl)*dz enddo open(3,file='lwpm') do k=1,kkt write(3,101) lwpm(k),zz(k) enddo close(3) open(3,file='lwpv') do k=1,kkt write(3,101) lwpv(k),zz(k) enddo close(3) open(3,file='iwpm') do k=1,kkt write(3,101) iwpm(k),zz(k) enddo close(3) open(3,file='iwpv') do k=1,kkt write(3,101) iwpv(k),zz(k) enddo close(3) c.... Calculate layer optical depth c... The approximation used: c tau = 2*pi*(3*rhoa/(pi*rhow))**(2/3) *h *Nc**(1/3) *(rc)**(2/3) c where Nc = particle concentration c... cloud water tNc = (2.e8)**0.3333333 pi=3.14159 rhow=1000. kk=0 do k=ks,ke,kl kk=kk+1 taucm(kk)=0.0 do l=k,k+kl do i=is,ie tau = 2.*pi*(3.*rho(k)/(pi*rhow))**(2./3.) *dz & *tNc*(qc(i,l)) **(2./3.) tauc(i,l)=tau taucm(kk)=tau+taucm(kk) enddo enddo taucm(kk)=taucm(kk)/tot enddo kk=0 do k=ks,ke,kl kk=kk+1 taucv(kk)=0.0 do l=k,k+kl do i=is,ie tau = 2.*pi*(3.*rho(k)/(pi*rhow))**(2./3.) *dz & *tNc*(qc(i,l)) **(2./3.) taucv(kk)=(tau-taucm(kk))**2 + taucv(kk) enddo enddo taucv(kk)=taucv(kk)/tot enddo open(3,file='taucm') do k=1,kkt write(3,101) taucm(k),zz(k) enddo close(3) open(3,file='taucv') do k=1,kkt write(3,101) taucv(k),zz(k) enddo close(3) c.... ice pi=3.14159 rhow=900. kk=0 do k=ks,ke,kl kk=kk+1 tauim(kk)=0.0 do l=k,k+kl do i=is,ie tNc=(rho(k)*qi(i,l)/5.e-9)**0.33333333 tau = 2.*pi*(3.*rho(k)/(pi*rhow))**(2./3.) *dz & *tNc*(qc(i,l)) **(2./3.) c tau = 2*pi*(3*rho(k)/(pi*rhow))**(2/3) *dz *tNc*(qi(i,l)**(2/3) taui(i,l)=tau tauim(kk)=tau+tauim(kk) enddo enddo taucm(kk)=taucm(kk)/tot enddo kk=0 do k=ks,ke,kl kk=kk+1 tauiv(kk)=0.0 do l=k,k+kl do i=is,ie tNc=(rho(k)*qi(i,l)/5.e-9)**0.33333333 tau = 2.*pi*(3.*rho(k)/(pi*rhow))**(2./3.) *dz & *tNc*(qi(i,l)) **(2./3.) c tau = 2*pi*(3*rho(k)/(pi*rhow))**(2/3) *dz *tNc*(qi(i,l)**(2/3) tauiv(kk)=(tau-tauim(kk))**2 + tauiv(kk) enddo enddo tauiv(kk)=tauiv(kk)/tot enddo open(3,file='tauim') do k=1,kkt write(3,101) tauim(k),zz(k) enddo close(3) open(3,file='tauiv') do k=1,kkt write(3,101) tauiv(k),zz(k) enddo close(3) c.... Cloud coverage ... fkl=float(kl+1) kk=0 ds = 1./float(ie-is+1) do k=ks,ke,kl kk=kk+1 sigma(kk)=0. do i=is,ie qcmm=0. qimm=0. do l=k,k+kl qcmm=qcmm+qc(i,l)/fkl qimm=qimm+qi(i,l)/fkl enddo if(qcmm+qimm.gt.1.e-6) then sigma(kk)=sigma(kk)+ds endif enddo enddo open(3,file='totcldf') do k=1,kkt write(3,101) sigma(k),zz(k) enddo close(3) c.... Convective cloud fraction .. kk=0 ds = 1./float(ie-is+1) do k=ks,ke,kl kk=kk+1 sigmac(kk)=0. do i=is,ie qcmm=0. qimm=0. wmmm=0. do l=k,k+kl qcmm=qcmm+qc(i,l)/fkl qimm=qimm+qi(i,1)/fkl wmmm=w(i,l)+w(i,l)/fkl enddo if(qcmm+qimm.gt.2.e-4 .and. wmmm.ge. 0.4) then icloud(i)=1 else icloud(i)=0 endif enddo do i=is,ie if(icloud(i).eq.1) then sigmac(kk)=sigmac(kk)+ds endif enddo enddo open(3,file='concldf') do k=1,kkt write(3,101) sigmac(k),zz(k) enddo close(3) c... Calculation of cloud coverage for ISCCP comparison c... calculate column optical depth do i=is,ie taul(i)=0. enddo do k=ks,ke do i=is,ie taul(i)=taui(i,k)+tauc(i,k) enddo enddo ds = 1./float(ie-is+1) sigmal=0. do i=is,ie if(taul(i).ge.0.1) then sigmal = sigmal+ds endif enddo write(2,203) sigmal 203 format('cloud fraction using optical depth threshold= ',e16.9) c.... Calculate eddy fluxes of moenentum, heat and moisture ds = 1./float(ie-is+1) C... mean u do k=ks,ke um(k)=0. do i=is,ie um(k)=um(k)+u(i,k) enddo um(k)=um(k)*ds enddo open(3,file='umean') do k=ks,ke write(3,101) um(k),z(k) enddo close(3) C... mean v do k=ks,ke vm(k)=0. do i=is,ie vm(k)=vm(k)+v(i,k) enddo vm(k)=vm(k)*ds enddo open(3,file='vmean') do k=ks,ke write(3,101) vm(k),z(k) enddo close(3) C... mean w do k=ks,ke wm(k)=0. do i=is,ie wm(k)=wm(k)+w(i,k) enddo wm(k)=wm(k)*ds enddo open(3,file='wmean') do k=ks,ke write(3,101) wm(k),z(k) enddo close(3) C... mean Theta do k=ks,ke thm(k)=0. do i=is,ie thm(k)=thm(k)+th(i,k) enddo thm(k)=thm(k)*ds enddo C... mean water vapor mixing ratio do k=ks,ke rvm(k)=0. do i=is,ie rvm(k)=rvm(k)+qv(i,k) enddo rvm(k)=rvm(k)*ds enddo c... Calculate eddy fluxes c... vertical u fluxes do k=ks,ke ueddy(k)=0.0 do i=is,ie ueddy(k)=(w(i,k)-wm(k))*(u(i,k)-um(k))+ueddy(k) enddo ueddy(k)=rho(k)*ueddy(k)*ds enddo open(3,file='ueddyf') do k=ks,ke write(3,101) ueddy(k),z(k) enddo close(3) c... v fluxes do k=ks,ke veddy(k)=0.0 do i=is,ie veddy(k)=(w(i,k)-vm(k))*(v(i,k)-vm(k))+veddy(k) enddo veddy(k)=rho(k)*veddy(k)*ds enddo open(3,file='veddyf') do k=ks,ke write(3,101) veddy(k),z(k) enddo close(3) c... w fluxes do k=ks,ke weddy(k)=0.0 do i=is,ie weddy(k)=(w(i,k)-wm(k))*(w(i,k)-wm(k))+weddy(k) enddo weddy(k)=rho(k)*weddy(k)*ds enddo open(3,file='weddyf') do k=ks,ke write(3,101) weddy(k),z(k) enddo close(3) c... theta fluxes do k=ks,ke theddy(k)=0.0 do i=is,ie theddy(k)=(w(i,k)-wm(k))*(th(i,k)-thm(k))+theddy(k) enddo theddy(k)=rho(k)*theddy(k)*ds enddo open(3,file='theddyf') do k=ks,ke write(3,101) theddy(k),z(k) enddo close(3) c... moisture fluxes do k=ks,ke rveddy(k)=0.0 do i=is,ie rveddy(k)=(w(i,k)-wm(k))*(qv(i,k)-qvm(k))+rveddy(k) enddo rveddy(k)=rho(k)*rveddy(k)*ds enddo open(3,file='vapeddyf') do k=ks,ke write(3,101) rveddy(k),z(k) enddo close(3) c.... Calculate apparent heat source Q1, neglecting diffusion ds=float(ie-is+1) print*,'ds = ',ds print*,'ks+1,ke-1= ',ks+1,ke-1 print*,'is,ie= ',is,ie do k=ks+1,ke-1 wdd=0. qd = 0. do i=is,ie d1=rho(k-1)*(w(i,k-1)-wm(k-1))*(th(i,k-1)-thm(k-1)) d2=rho(k+1)*(w(i,k+1)-wm(k+1))*(th(i,k+1)-thm(k+1)) wdd = wdd + (d2-d1)*0.5/dz qd = qd + qdot(i,k) enddo wdd=wdd/ds qd=qd/ds q1(k) = -(1/(pe(k)*rho(k))) * wdd + qd print*,'k= ',k,' qdot =', qd*3600. sw1(k) = qd sw2(k) = q1(k) - qd enddo open(3,file='q1') do k=ks+1,ke-1 write(3,101) 3600.*q1(k),z(k) enddo close(3) open(3,file='qdd') do k=ks+1,ke-1 write(3,101) 3600.*sw1(k),z(k) enddo close(3) open(3,file='qwd') do k=ks+1,ke-1 write(3,101) 3600.*sw2(k),z(k) enddo close(3) c.... Calculate apparent moisture sink Q2, neglecting diffusion cons=2.5e6/1005.0 ds=float(ie-is+1) do k=ks+1,ke-1 wdd=0. qd = 0. do i=is,ie d1=rho(k-1)*(w(i,k-1)-wm(k-1))*(qv(i,k-1)-rvm(k-1)) d2=rho(k+1)*(w(i,k+1)-wm(k+1))*(qv(i,k+1)-rvm(k+1)) wdd = wdd + (d2-d1)*0.5/dz qd = qd + qdot2(i,k) enddo wdd=wdd/ds qd=qd/ds q2(k) = (cons/rho(k)) * wdd + qd enddo open(3,file='q2') do k=ks+1,ke-1 write(3,101) q2(k),z(k) enddo close(3) c.... Calculate apparent u momentum source Qu, neglecting diffusion ds=float(ie-is+1) do k=ks+1,ke-1 wdd=0. qd = 0. do i=is,ie d1=rho(k-1)*(w(i,k-1)-wm(k-1))*(u(i,k-1)-um(k-1)) d2=rho(k+1)*(w(i,k+1)-wm(k+1))*(u(i,k+1)-um(k+1)) wdd = wdd + (d2-d1)*0.5/dz enddo wdd=wdd/ds qu(k) = -(1.0/rho(k)) * wdd enddo open(3,file='qu') do k=ks+1,ke-1 write(3,101) qu(k),zz(k) enddo close(3) c.... Calculate apparent v momentum source Qv, neglecting diffusion ds=float(ie-is+1) do k=ks+1,ke-1 wdd=0. qd = 0. do i=is,ie d1=rho(k-1)*(w(i,k-1)-wm(k-1))*(v(i,k-1)-vm(k-1)) d2=rho(k+1)*(w(i,k+1)-wm(k+1))*(v(i,k+1)-vm(k+1)) wdd = wdd + (d2-d1)*0.5/dz enddo wdd=wdd/ds qqv(k) = -(1.0/rho(k)) * wdd enddo open(3,file='qv') do k=ks+1,ke-1 write(3,101) qqv(k),zz(k) enddo close(3) stop end