Blame | Last modification | View Log | Download | RSS feed
subroutine zlev(z,t,q,zb,sp,ie,je,ke,ak,bk)c ===========================================c argument declarationinteger ie,je,ke,its,iqsreal z(ie,je,ke),t(ie,je,ke),getpsreal q(ie,je,ke),zb(ie,je),sp(ie,je)real ak(ke),bk(ke)c variable declarationinteger i,j,kreal pu, po, zo, zu, tvreal rdg,tzero,epsdata rdg,tzero,eps /29.271, 273.15, 0.6078/c t and q are staggered (stag(3)=-0.5; its=int(2*stag(3))its=-1iqs=-1c computation of height of main levelsdo i=1,iedo j=1,jezu = zb(i,j)pu = sp(i,j)tv = (t(i,j,1)+tzero)*(1.+eps*q(i,j,1))po = ak(1)+bk(1)*sp(i,j)if (po.le.0) po=0.5*puzo = zu + rdg*tv*alog(pu/po)z(i,j,1) = zopu = pozu = zodo k=2,ketv = (0.5*(t(i,j,k)+t(i,j,k-1-its))+tzero)*> (1.+eps*0.5*(q(i,j,k)+q(i,j,k-1-iqs)))po = ak(k)+bk(k)*sp(i,j)if (po.le.0) po=0.5*puzo = zu + rdg*tv*alog(pu/po)z(i,j,k) = zopu = pozu = zoenddoenddoenddoendsubroutine flightlev(fl,p,ie,je,ke)c ===================================c Convert an array of pressures (hPa) into an array of altitudes (m)c or vice versa using the ICAO standard atmosphere.cc See http://www.pdas.com/coesa.htmcc The 7 layers of the US standard atmosphere are:cc h1 h2 dT/dh h1,h2 geopotential alt in kmc 0 11 -6.5 dT/dh in K/kmc 11 20 0.0c 20 32 1.0c 32 47 2.8c 47 51 0.0c 51 71 -2.8c 71 84.852 -2.0c based upon IDL routine from Dominik Brunnerc converted to f90: Jan 2002 Heini Wernliimplicit nonec argument declarationinteger ie,je,kereal fl(ie,je,ke),p(ie,je,ke)c variable declarationinteger i,j,k,n,ilreal limits(0:7),lrs(7)integer iszero(7)real pb(0:7),tb(0:7)real g,r,gmr,feetdata g,r,feet /9.80665, 287.053, 0.3048/c define layer boundaries in m, lapsre rates (9 means 0) and isothermal flaglimits(0)=0.limits(1)=11000.limits(2)=20000.limits(3)=32000.limits(4)=47000.limits(5)=51000.limits(6)=71000.limits(7)=84852.lrs(0)=-6.5/1000.lrs(1)=9.0/1000.lrs(2)=1.0/1000.lrs(3)=2.8/1000.lrs(4)=9.0/1000.lrs(5)=-2.8/1000.lrs(6)=-2.0/1000.iszero(0)=0iszero(1)=1iszero(2)=0iszero(3)=0iszero(4)=1iszero(5)=0iszero(6)=0gmr=g/r ! Hydrostatic constantc calculate pressures at layer boundariespB(0)=1013.25 ! pressure at surfaceTB(0)=288.15 ! Temperature at surfacec loop over layers and get pressures and temperatures at level topsdo i=0,6TB(i+1)=TB(i)+(1-iszero(i))*lrs(i)*(limits(i+1)-limits(i))pB(i+1)=(1-iszero(i))*pB(i)*exp(alog(TB(i)/TB(i+1))*gmr/lrs(i))+> iszero(i)*PB(i)*exp(-gmr*(limits(i+1)-limits(i))/TB(i))enddoc now calculate which layer each value belongs toc and calculate the corresponding altitudesdo i=1,iedo j=1,jedo k=1,kedo n=0,6if ((pb(n).ge.p(i,j,k)).and.(pb(n+1).le.p(i,j,k))) thenil=ngoto 100endifenddo100 continuefl(i,j,k)=iszero(il)*(-ALOG(p(i,j,k)/pB(il))*TB(il)/gmr+> limits(il))+(1-iszero(il))*((TB(il)/((p(i,j,k)/> pB(il))**(lrs(il)/gmr))-TB(il))/lrs(il)+limits(il))fl(i,j,k)=fl(i,j,k)/feet/100.enddoenddoenddoendsubroutine pottemp(pt,t,sp,ie,je,ke,ak,bk)c ==========================================c argument declarationinteger ie,je,kereal pt(ie,je,ke),t(ie,je,ke),sp(ie,je),> ak(ke),bk(ke)c variable declarationinteger i,j,kreal rdcp,tzerodata rdcp,tzero /0.286,273.15/c statement-functions for the computation of pressurereal pp,psrfinteger ispp(is)=ak(is)+bk(is)*psrfc computation of potential temperaturedo i=1,iedo j=1,jepsrf=sp(i,j)do k=1,kec distinction of temperature in K and deg Cif (t(i,j,k).lt.100.) thenpt(i,j,k)=(t(i,j,k)+tzero)*( (1000./pp(k))**rdcp )elsept(i,j,k)=t(i,j,k)*( (1000./pp(k))**rdcp )endifenddoenddoenddoendsubroutine tnat(tn,t,p,ie,je,ke)c ================================c argument declarationinteger ie,je,kereal tn(ie,je,ke),t(ie,je,ke),p(ie,je,ke)c variable declarationinteger i,j,kreal tzerodata tzero /273.15/real mb_to_mmreal a1,a2,a3,a4,a5,c0,c1,c2,dreal chiH2O,chiHNO3,pH2O,pHNO3,log_H2O,log_NATmb_to_mm=760./1000.a1 = -2.7836a2 = -0.00088a3 = 89.7674a4 = -26242.0a5 = 0.021135chiH2O = 5.e-6chiHNO3 = 5.e-9do i=1,iedo j=1,jedo k=1,kepH2O = p(i,j,k)*chiH2OpHNO3 = p(i,j,k)*chiHNO3log_H2O = log(pH2O*mb_to_mm)log_NAT = log(pHNO3*mb_to_mm)c0 = log_NAT-log_H2O*a1-a3c1 = a2*log_H2O+a5c2 = a4d = c0*c0-4.0*c1*c2tn(i,j,k)=t(i,j,k)+tzero-(c0+sqrt(d))/(2.*c1)enddoenddoenddoendsubroutine read_tice(tice_arr)c ==============================integer nreal tice_arr(100)open(10,>file='/home/henry/prog/tnat_tice/t_ice_from_10_every_2_hPa')do n=1,71read(10,*)tice_arr(n)enddoclose(10)endsubroutine tice(ti,t,p,tice_arr,ie,je,ke)c =========================================c argument declarationinteger ie,je,kereal ti(ie,je,ke),t(ie,je,ke),p(ie,je,ke)c variable declarationinteger i,j,kreal tzerodata tzero /273.15/real tice_arr(100)real p0,p1,pinc,rindinteger ind,npc define the lowest p-value and the p-increment for the tablep0=10.pinc=2.np=71p1=p0+real(np-1)*pincdo i=1,iedo j=1,jedo k=1,keif (p(i,j,k).lt.p0) thenprint*,'*** error: table is not prepared for p < p0 ***'print*,ie,je,ke,i,j,k,p(i,j,k),p0call exit(1)endifif (p(i,j,k).gt.p1) thenprint*,'*** error: table is not prepared for p > p1 ***'call exit(1)endifc calculate the real index of the given pressure levelrind=(p(i,j,k)-p0)/pinc+1.ind=int(rind)rind=rind-real(ind)c interpolate ticeti(i,j,k)=t(i,j,k)+tzero-> (tice_arr(ind)+rind*(tice_arr(ind+1)-tice_arr(ind)))enddoenddoenddoendsubroutine temp(t,pt,p,ie,je,ke,ak,bk)c ======================================c argument declarationinteger ie,je,kereal pt(ie,je,ke),t(ie,je,ke),p(ie,je,ke),> ak(ke),bk(ke)c variable declarationinteger i,j,kreal rdcp,tzerodata rdcp,tzero /0.286,273.15/c computation of potential temperaturedo i=1,iedo j=1,jedo k=1,ket(i,j,k)=pt(i,j,k)*( (p(i,j,k)/1000.)**rdcp )enddoenddoenddoendsubroutine laitpv(lpv,pv,th,ie,je,ke)c =====================================c argument declarationinteger ie,je,kereal lpv(ie,je,ke),pv(ie,je,ke),th(ie,je,ke)c variable declarationinteger i,j,kreal rdcp,tzerodata rdcp,tzero /0.286,273.15/c computation of Lait PVdo i=1,iedo j=1,jedo k=1,kelpv(i,j,k)=pv(i,j,k)*((th(i,j,k)/420.)**(-9./2.))enddoenddoenddoendsubroutine relhum(rh,q,t,sp,ie,je,ke,ak,bk)c ===========================================c argument declarationinteger ie,je,kereal rh(ie,je,ke),t(ie,je,ke),q(ie,je,ke),> sp(ie,je),ak(ke),bk(ke)c variable declarationinteger i,j,kreal rdcp,tzeroreal b1,b2w,b3,b4w,r,rd,gqd,gedata rdcp,tzero /0.286,273.15/data b1,b2w,b3,b4w,r,rd /6.1078, 17.2693882, 273.16, 35.86,& 287.05, 461.51/c statement-functions for the computation of pressurereal pp,psrfinteger ispp(is)=ak(is)+bk(is)*psrfdo i=1,iedo j=1,jepsrf=sp(i,j)do k=1,kege = b1*exp(b2w*t(i,j,k)/(t(i,j,k)+b3-b4w))gqd= r/rd*ge/(pp(k)-(1.-r/rd)*ge)rh(i,j,k)=100.*q(i,j,k)/gqdenddoenddoenddoendsubroutine relhum_i(rh,q,p,ie,je,ke,ak,mdv)c ===========================================c argument declarationinteger ie,je,kereal rh(ie,je,ke),p(ie,je,ke),q(ie,je,ke),> ak(ke),mdvc variable declarationinteger i,j,kreal p0,kappadata p0,kappa /1000.,0.286/real rdcp,tzerodata rdcp,tzero /0.286,273.15/real b1,b2w,b3,b4w,r,rd,gqd,gedata b1,b2w,b3,b4w,r,rd /6.1078, 17.2693882, 273.16, 35.86,& 287.05, 461.51/do i=1,iedo j=1,jedo k=1,keif (p(i,j,k).eq.mdv) thenrh(i,j,k)=mdvelsett=ak(k)*((p(i,j,k)/p0)**kappa)-tzeroge = b1*exp(b2w*tt/(tt+b3-b4w))gqd= r/rd*ge/(p(i,j,k)-(1.-r/rd)*ge)rh(i,j,k)=100.*q(i,j,k)/gqdendifenddoenddoenddoendsubroutine gradth(gth,th,sp,levtyp,cl,ie,je,ke,ak,bk,vmin,vmax)C ===============================================================c argument declarationinteger ie,je,ke,levtypreal gth(ie,je,ke),th(ie,je,ke),sp(ie,je),cl(ie,je)real ak(ke),bk(ke),vmin(4),vmax(4)c variable declarationREAL,ALLOCATABLE, DIMENSION (:,:) :: dspdx,dspdyREAL,ALLOCATABLE, DIMENSION (:,:,:) :: dthdx,dthdyinteger i,j,k,ind,ind2,statallocate(dthdx(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'gradth: error allocating dthdx'allocate(dthdy(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'gradth: error allocating dthdy'allocate(dspdx(ie,je),STAT=stat)if (stat.ne.0) print*,'gradth: error allocating dspdx'allocate(dspdy(ie,je),STAT=stat)if (stat.ne.0) print*,'gradth: error allocating dspdy'if (levtyp.eq.1) thencall ddh2(th,dthdx,cl,'X',ie,je,1,vmin,vmax)call ddh2(th,dthdy,cl,'Y',ie,je,1,vmin,vmax)else if (levtyp.eq.3) thencall ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)call ddh3(th,dthdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)call ddh3(th,dthdy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)endifgth=sqrt(dthdx**2.+dthdy**2.)IF (ALLOCATED(dspdx)) DEALLOCATE(dspdx)IF (ALLOCATED(dspdy)) DEALLOCATE(dspdy)IF (ALLOCATED(dthdx)) DEALLOCATE(dthdx)IF (ALLOCATED(dthdy)) DEALLOCATE(dthdy)endsubroutine calc_gradpv(gpv,pv,cl,ie,je,ke,ak,bk,vmin,vmax)C ===============================================================c argument declarationinteger ie,je,ke,levtypreal gpv(ie,je,ke),pv(ie,je,ke),sp(ie,je),cl(ie,je)real ak(ke),bk(ke),vmin(4),vmax(4)c variable declarationREAL,ALLOCATABLE, DIMENSION (:,:,:) :: dpvdx,dpvdyinteger i,j,k,ind,ind2,statallocate(dpvdx(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'gradpv: error allocating dpvdx'allocate(dpvdy(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'gradpv: error allocating dpvdy'call ddh2(pv,dpvdx,cl,'X',ie,je,1,vmin,vmax)call ddh2(pv,dpvdy,cl,'Y',ie,je,1,vmin,vmax)gpv=1e6*sqrt(dpvdx**2.+dpvdy**2.)IF (ALLOCATED(dpvdx)) DEALLOCATE(dpvdx)IF (ALLOCATED(dpvdy)) DEALLOCATE(dpvdy)endsubroutine wetbpt(thw,the,sp,ie,je,ke,ak,bk)c ============================================c argument declarationinteger ie,je,kereal thw(ie,je,ke),the(ie,je,ke),> sp(ie,je),ak(ke),bk(ke)c variable declarationreal tsainteger i,j,kdo k=1,kedo j=1,jedo i=1,iethw(i,j,k)=tsa(the(i,j,k)-273.15,1000.)+273.15enddoenddoenddoendreal function tsa(os,p)c =======================C This function returns the temperature tsa (celsius) on ac saturation adiabat at pressure p (millibars). os is the equivalentc potential temperature of the parcel (celsius). sign(a,b) replacesc the algebraic sign of a with that of b.C b is an empirical constant approximately equal to 0.001 of thec latent heat of vaporization for water divided by the specificc heat at constant pressure for dry air.real a,b,os,p,tq,d,tqk,x,winteger idata b/2.6518986/a=os+273.16C tq is the first guess for tsatq=253.16C d is an initial value used in the iteration belowd=120.C Iterate to obtain sufficient accuracy....see table 1, p.8C of Stipanuk (1973) for equation used in iterationdo 1 i=1,12tqk=tq-273.16d=d/2.x=a*exp(-b*w(tqk,p)/tq)-tq*((1000./p)**.286)if (abs(x).lt.1e-7) go to 2tq=tq+sign(d,x)1 continue2 tsa=tq-273.16endreal function w(t,p)c ====================C This function returns the mixing ratio (grams of water vapor perC kilogram of dry air) given the dew point (celsius) and pressureC (millibars). If the temperature is input instead of theC dew point, then saturation mixing ratio (same units) is returned.C The formula is found in most meteorological texts.real t,p,tkel,x,esattkel=t+273.16x=esat(tkel) ! our function esat requires t in Kelvinw=622.*x/(p-x)endreal function esat(t)c =====================C This function returns the saturation vapor pressure over waterc (mb) given the temperature (Kelvin).C The algorithm is due to Nordquist, W. S. ,1973: "NumericalC Approximations of Selected Meteorological Parameters for CloudC Physics Problems" ECOM-5475, Atmospheric Sciences Laboratory,c U. S. Army Electronics Command, White Sands Missile Range,c New Mexico 88002.real p1,p2,c1,tp1=11.344-0.0303998*tp2=3.49149-1302.8844/tc1=23.832241-5.02808*log10(t)esat=10.**(c1-1.3816e-7*10.**p1+8.1328e-3*10.**p2-2949.076/t)endsubroutine equpot(ap,t,q,sp,ie,je,ke,ak,bk)c ===========================================c calculate equivalent potential temperaturec argument declarationinteger ie,je,kereal ap(ie,je,ke),t(ie,je,ke),sp(ie,je)real q(ie,je,ke),ak(ke),bk(ke)c variable declarationinteger i,j,kreal rdcp,tzerodata rdcp,tzero /0.286,273.15/c statement-functions for the computation of pressurereal pp,psrfinteger ispp(is)=ak(is)+bk(is)*psrfc computation of potential temperaturedo i=1,iedo j=1,jepsrf=sp(i,j)do k=1,keap(i,j,k) = (t(i,j,k)+tzero)*(1000./pp(k))+ **(0.2854*(1.0-0.28*q(i,j,k)))*exp(+ (3.376/(2840.0/(3.5*alog(t(i,j,k)+tzero)-alog(+ 100.*pp(k)*max(1.0E-10,q(i,j,k))/(0.622+0.378*+ q(i,j,k)))-0.1998)+55.0)-0.00254)*1.0E3*+ max(1.0E-10,q(i,j,k))*(1.0+0.81*q(i,j,k)))enddoenddoenddoendsubroutine diabheat(dhr,t,w,rh,sp,ie,je,ke,ak,bk)c =================================================c argument declarationinteger ie,je,kereal dhr(ie,je,ke),t(ie,je,ke),w(ie,je,ke),& rh(ie,je,ke),sp(ie,je),ak(ke),bk(ke)c variable declarationinteger i,j,kreal p0,kappa,tzerodata p0,kappa,tzero /1000.,0.286,273.15/real blog10,cp,r,lw,epsdata blog10,cp,r,lw,eps /.08006,1004.,287.,2.5e+6,0.622/real esat,c,ttc statement-functions for the computation of pressurereal pp,psrfinteger ispp(is)=ak(is)+bk(is)*psrfc computation of diabatic heating ratedo i=1,iedo j=1,jepsrf=sp(i,j)do k=1,keif (rh(i,j,k).lt.80.) then ! only moist air of interestdhr(i,j,k)=0. ! cond. heating rate set to zeroelse if (w(i,j,k).gt.0.) then ! cond. heating only! for ascentdhr(i,j,k)=0.elsett=t(i,j,k)*((pp(k)/p0)**kappa) ! temp. from pot.temp.c=lw/cp*eps*blog10*esat(tt)/pp(k)dhr(i,j,k)=21600.* ! in units K per 6 hours> (1.-exp(.2*(80.-rh(i,j,k)))) ! weighting fun.! for 80<RH<100> *(-c*kappa*t(i,j,k)*w(i,j,k)/(100.*pp(k)))/(1.+c)endifenddoenddoenddoendsubroutine diabheat2(dhr,t,w,rh,sp,ie,je,ke,ak,bk)c ==================================================c argument declarationinteger ie,je,kereal dhr(ie,je,ke),t(ie,je,ke),w(ie,je,ke),& rh(ie,je,ke),sp(ie,je),ak(ke),bk(ke)c variable declarationinteger i,j,kreal p0,kappa,tzerodata p0,kappa,tzero /1000.,0.286,273.15/real blog10,cp,r,lw,epsdata blog10,cp,r,lw,eps /.08006,1004.,287.,2.5e+6,0.622/real esat,c,ttc statement-functions for the computation of pressurereal pp,psrfinteger ispp(is)=ak(is)+bk(is)*psrfc computation of diabatic heating ratedo i=1,iedo j=1,jepsrf=sp(i,j)do k=1,keif (rh(i,j,k).lt.80.) then ! only moist air of interestdhr(i,j,k)=0. ! cond. heating rate set to zeroelse if (w(i,j,k).gt.0.) then ! cond. heating only! for ascentdhr(i,j,k)=0.elsec=lw/cp*eps*blog10*esat((t(i,j,k)+273.15))/pp(k)tt=(t(i,j,k)+273.15)*((p0/pp(k))**kappa) ! pot.temp from tempdhr(i,j,k)=21600.* ! in units K per 6 hours> (1.-exp(.2*(80.-rh(i,j,k)))) ! weighting fun.! for 80<RH<100> *(-c*kappa*tt*w(i,j,k)/(100.*pp(k)))/(1.+c)endifenddoenddoenddoendsubroutine diabpvr(dpvr,uu,vv,dhr,sp,cl,f,ie,je,ke,ak,bk,> vmin,vmax)C =========================================================c argument declarationinteger ie,je,kereal,intent(OUT) :: dpvr(ie,je,ke)real,intent(IN) :: uu(ie,je,ke),vv(ie,je,ke),> dhr(ie,je,ke),sp(ie,je),cl(ie,je),f(ie,je)real,intent(IN) :: ak(ke),bk(ke),vmin(4),vmax(4)c variable declarationREAL,ALLOCATABLE, DIMENSION (:,:) :: dspdx,dspdyREAL,ALLOCATABLE, DIMENSION (:,:,:) :: dhrdp,dudp,dvdp,dvdxREAL,ALLOCATABLE, DIMENSION (:,:,:) :: dhrdx,dudy,dhrdyinteger k,statallocate(dspdx(ie,je),STAT=stat)if (stat.ne.0) print*,'diabpvr: error allocating dspdx'allocate(dspdy(ie,je),STAT=stat)if (stat.ne.0) print*,'diabpvr: error allocating dspdy'allocate(dhrdp(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'diabpvr: error allocating dhrdp'allocate(dudp(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'diabpvr: error allocating dudp'allocate(dvdp(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'diabpvr: error allocating dvdp'allocate(dvdx(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'diabpvr: error allocating dvdx'allocate(dhrdx(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'diabpvr: error allocating dhrdx'allocate(dudy(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'diabpvr: error allocating dudy'allocate(dhrdy(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'diabpvr: error allocating dhrdy'call ddp(dhr,dhrdp,sp,ie,je,ke,ak,bk)call ddp(uu,dudp,sp,ie,je,ke,ak,bk)call ddp(vv,dvdp,sp,ie,je,ke,ak,bk)call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)call ddh3(dhr,dhrdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)call ddh3(vv,dvdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)call ddh3(dhr,dhrdy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)call ddh3(uu,dudy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)do k=1,kedpvr(1:ie,1:je,k)=-1.E6*9.80665*(> -dvdp(1:ie,1:je,k)*dhrdx(1:ie,1:je,k)> +dudp(1:ie,1:je,k)*dhrdy(1:ie,1:je,k)> +(-dudy(1:ie,1:je,k)+dvdx(1:ie,1:je,k)> +f(1:ie,1:je))*dhrdp(1:ie,1:je,k))enddoIF (ALLOCATED(dspdx)) DEALLOCATE(dspdx)IF (ALLOCATED(dspdy)) DEALLOCATE(dspdy)IF (ALLOCATED(dhrdp)) DEALLOCATE(dhrdp)IF (ALLOCATED(dudp)) DEALLOCATE(dudp)IF (ALLOCATED(dvdp)) DEALLOCATE(dvdp)IF (ALLOCATED(dvdx)) DEALLOCATE(dvdx)IF (ALLOCATED(dhrdx)) DEALLOCATE(dhrdx)IF (ALLOCATED(dudy)) DEALLOCATE(dudy)IF (ALLOCATED(dhrdy)) DEALLOCATE(dhrdy)endsubroutine vort(var,uu,vv,sp,cl,ie,je,ke,ak,bk,vmin,vmax)C =========================================================C calculate vorticity VORT=D[V:X]-D[U*COS:Y]/COSc argument declarationinteger ie,je,kereal,intent(OUT) :: var(ie,je,ke)real,intent(IN) :: uu(ie,je,ke),vv(ie,je,ke),> sp(ie,je),cl(ie,je)real ak(ke),bk(ke),vmin(4),vmax(4)c variable declarationREAL,ALLOCATABLE, DIMENSION (:,:) :: dspdx,dspdyREAL,ALLOCATABLE, DIMENSION (:,:,:) :: dvdx,dudy,uuclinteger statreal poleallocate(dspdx(ie,je),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dspdx(ie,je)'allocate(dspdy(ie,je),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dspdy(ie,je)'allocate(dvdx(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dvdx'allocate(uucl(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'vort: error allocating uucl'allocate(dudy(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dudy'call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)call ddh3(vv,dvdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)do k=1,keuucl(1:ie,1:je,k)=uu(1:ie,1:je,k)*cl(1:ie,1:je)enddocall ddh3(uucl,dudy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)do k=1,kevar(1:ie,1:je,k)=1.E4*(dvdx(1:ie,1:je,k)-> dudy(1:ie,1:je,k)/cl(1:ie,1:je))* var(1:ie,1:je,k)=1.E4*(dvdx(1:ie,1:je,k)-* > dudy(1:ie,1:je,k))enddoc interpolate poles from neighbouring points on every levelif (vmax(2).eq.90.) thendo k=1,kepole=0.do i=1,iepole=pole+var(i,je-1,k)enddopole=pole/real(ie)do i=1,ievar(i,je,k)=poleenddoenddoendifif (vmin(2).eq.-90.) thendo k=1,kepole=0.do i=1,iepole=pole+var(i,2,k)enddopole=pole/real(ie)do i=1,ievar(i,1,k)=poleenddoenddoendifIF (ALLOCATED(dspdx)) DEALLOCATE(dspdx)IF (ALLOCATED(dspdy)) DEALLOCATE(dspdy)IF (ALLOCATED(dvdx)) DEALLOCATE(dvdx)IF (ALLOCATED(dudy)) DEALLOCATE(dudy)IF (ALLOCATED(uucl)) DEALLOCATE(uucl)endsubroutine avort(var,uu,vv,sp,cl,f,ie,je,ke,ak,bk,vmin,vmax)C ============================================================C calculate absolute vorticity AVO=F+D[V:X]-D[U*COS:Y]/COSc argument declarationinteger ie,je,kereal,intent(OUT) :: var(ie,je,ke)real,intent(IN) :: uu(ie,je,ke),vv(ie,je,ke),> sp(ie,je),cl(ie,je),f(ie,je)real ak(ke),bk(ke),vmin(4),vmax(4)c variable declarationREAL,ALLOCATABLE, DIMENSION (:,:) :: dspdx,dspdyREAL,ALLOCATABLE, DIMENSION (:,:,:) :: dvdx,dudy,uuclinteger statreal poleallocate(dspdx(ie,je),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dspdx(ie,je)'allocate(dspdy(ie,je),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dspdy(ie,je)'allocate(dvdx(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dvdx'allocate(uucl(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'vort: error allocating uucl'allocate(dudy(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dudy'call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)call ddh3(vv,dvdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)do k=1,keuucl(1:ie,1:je,k)=uu(1:ie,1:je,k)*cl(1:ie,1:je)enddocall ddh3(uucl,dudy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)do k=1,kevar(1:ie,1:je,k)=1.E4*(f(1:ie,1:je)+dvdx(1:ie,1:je,k)-> dudy(1:ie,1:je,k)/cl(1:ie,1:je))enddoc interpolate poles from neighbouring points on every levelif (vmax(2).eq.90.) thendo k=1,kepole=0.do i=1,iepole=pole+var(i,je-1,k)enddopole=pole/real(ie)do i=1,ievar(i,je,k)=poleenddoenddoendifif (vmin(2).eq.-90.) thendo k=1,kepole=0.do i=1,iepole=pole+var(i,2,k)enddopole=pole/real(ie)do i=1,ievar(i,1,k)=poleenddoenddoendifIF (ALLOCATED(dspdx)) DEALLOCATE(dspdx)IF (ALLOCATED(dspdy)) DEALLOCATE(dspdy)IF (ALLOCATED(dvdx)) DEALLOCATE(dvdx)IF (ALLOCATED(dudy)) DEALLOCATE(dudy)IF (ALLOCATED(uucl)) DEALLOCATE(uucl)endsubroutine vort_i(var,uu,vv,cl,ie,je,ke,vmin,vmax,mdv)C ======================================================C calculate vorticity VORT=1e4*D[V:X]-D[U*COS:Y]/COS on isentropicC levelsC with misdat (mdv) treatment !C mark liniger 990501c argument declarationinteger ie,je,kereal,intent(OUT) :: var(ie,je,ke)real,intent(IN) :: uu(ie,je,ke),vv(ie,je,ke),> cl(ie,je)real vmin(4),vmax(4),mdvc variable declarationREAL,ALLOCATABLE, DIMENSION (:,:) :: dvdx,dudy,uuclinteger statallocate(dvdx(ie,je),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dvdx'allocate(uucl(ie,je),STAT=stat)if (stat.ne.0) print*,'vort: error allocating uucl'allocate(dudy(ie,je),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dudy'do k=1,kecall ddh2m(vv(1,1,k),dvdx,cl,'X',ie,je,1,vmin,vmax,mdv)do i=1,iedo j=1,jeif (uu(i,j,k).eq.mdv) thenuucl(i,j)=mdvelseuucl(i,j)=uu(i,j,k)*cl(i,j)endifenddoenddocall ddh2m(uucl,dudy,cl,'Y',ie,je,1,vmin,vmax,mdv)do i=1,iedo j=1,jeif ((dvdx(i,j).eq.mdv).or.> (dudy(i,j).eq.mdv).or.> (cl(i,j).lt.1e-3)) thenvar(i,j,k)=mdvelsevar(i,j,k)=1.E4*(dvdx(i,j)-dudy(i,j)> /cl(i,j))endifenddoenddoenddoIF (ALLOCATED(dvdx)) DEALLOCATE(dvdx)IF (ALLOCATED(dudy)) DEALLOCATE(dudy)IF (ALLOCATED(uucl)) DEALLOCATE(uucl)endsubroutine curvo(var,uu,vv,sp,cl,ie,je,ke,ak,bk,vmin,vmax)C ==========================================================C calculate curvature vorticityC =u^2*d[v:x]-v^2*d[u*cos:y]/cos-u*v*(d[u:x]-d[v*cos:y]/cos))/(u^2+v^2)c argument declarationinteger ie,je,kereal,intent(OUT) :: var(ie,je,ke)real,intent(IN) :: uu(ie,je,ke),vv(ie,je,ke),> sp(ie,je),cl(ie,je)real ak(ke),bk(ke),vmin(4),vmax(4)c variable declarationREAL,ALLOCATABLE, DIMENSION (:,:) :: dspdx,dspdyREAL,ALLOCATABLE, DIMENSION (:,:,:) :: dudx,dvdx,dudy,dvdy,> uucl,vvclinteger statallocate(dspdx(ie,je),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dspdx'allocate(dspdy(ie,je),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dspdy'allocate(dudx(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dudx'allocate(dvdx(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dvdx'allocate(uucl(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'vort: error allocating uucl'allocate(vvcl(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'vort: error allocating vvcl'allocate(dudy(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dudy'allocate(dvdy(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'vort: error allocating dvdy'call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)call ddh3(uu,dudx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)call ddh3(vv,dvdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)do k=1,keuucl(1:ie,1:je,k)=uu(1:ie,1:je,k)*cl(1:ie,1:je)vvcl(1:ie,1:je,k)=vv(1:ie,1:je,k)*cl(1:ie,1:je)enddocall ddh3(uucl,dudy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)call ddh3(vvcl,dvdy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)do k=1,kedo j=1,jedo i=1,ievar(i,j,k)=1.e4*(uu(i,j,k)**2.*dvdx(i,j,k)-> vv(i,j,k)**2.*dudy(i,j,k)/cl(i,j)-> uu(i,j,k)*vv(i,j,k)*(dudx(i,j,k)-dvdy(i,j,k)/cl(i,j)))/> max(uu(i,j,k)**2.+vv(i,j,k)**2.,1.e-3)enddoenddoenddoc set poles values to mdvif (vmax(2).eq.90.) thendo k=1,kedo i=1,ievar(i,je,k)=-999.98999enddoenddoendifif (vmin(2).eq.-90.) thendo k=1,kedo i=1,ievar(i,1,k)=-999.98999enddoenddoendifIF (ALLOCATED(dspdx)) DEALLOCATE(dspdx)IF (ALLOCATED(dspdy)) DEALLOCATE(dspdy)IF (ALLOCATED(dudx)) DEALLOCATE(dudx)IF (ALLOCATED(dvdx)) DEALLOCATE(dvdx)IF (ALLOCATED(dudy)) DEALLOCATE(dudy)IF (ALLOCATED(dvdy)) DEALLOCATE(dvdy)IF (ALLOCATED(uucl)) DEALLOCATE(uucl)IF (ALLOCATED(vvcl)) DEALLOCATE(vvcl)endsubroutine potvort(pv,uu,vv,th,sp,cl,f,ie,je,ke,ak,bk,> vmin,vmax)C ======================================================C calculate PVc argument declarationinteger ie,je,kereal pv(ie,je,ke),uu(ie,je,ke),vv(ie,je,ke),> th(ie,je,ke),sp(ie,je),> cl(ie,je),f(ie,je)real ak(ke),bk(ke),vmin(4),vmax(4)real pvpolec variable declarationREAL,ALLOCATABLE, DIMENSION (:,:) :: dspdx,dspdyREAL,ALLOCATABLE, DIMENSION (:,:,:) :: dthdp,dudp,dvdpREAL,ALLOCATABLE, DIMENSION (:,:,:) :: uu2,dvdx,dudy,dthdx,dthdyinteger i,j,k,statallocate(dspdx(ie,je),STAT=stat)if (stat.ne.0) print*,'potvort: error allocating dspdx(ie,je)'allocate(dspdy(ie,je),STAT=stat)if (stat.ne.0) print*,'potvort: error allocating dspdy(ie,je)'allocate(dthdp(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort: error allocating dthdp'allocate(dudp(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort: error allocating dudp'allocate(dvdp(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort: error allocating dvdp'allocate(dvdx(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort: error allocating dvdx'allocate(dudy(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort: error allocating dudy'allocate(dthdx(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort: error allocating dthdx'allocate(dthdy(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort: error allocating dthdy'allocate(uu2(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort: error allocating uu2'call ddp(th,dthdp,sp,ie,je,ke,ak,bk)call ddp(uu,dudp,sp,ie,je,ke,ak,bk)call ddp(vv,dvdp,sp,ie,je,ke,ak,bk)call ddh2(sp,dspdx,cl,'X',ie,je,1,vmin,vmax)call ddh2(sp,dspdy,cl,'Y',ie,je,1,vmin,vmax)call ddh3(th,dthdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)call ddh3(vv,dvdx,sp,dspdx,cl,'X',ie,je,ke,vmin,vmax,ak,bk)call ddh3(th,dthdy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)c conversion of uu for spheric coordinates (uu*cos(phi))do k=1,keuu2(1:ie,1:je,k)=uu(1:ie,1:je,k)*cl(1:ie,1:je)enddocall ddh3(uu2,dudy,sp,dspdy,cl,'Y',ie,je,ke,vmin,vmax,ak,bk)c conversion of dudy for spheric coordinates (dudy/cos(phi))do k=1,kedudy(1:ie,1:je,k)=dudy(1:ie,1:je,k)/cl(1:ie,1:je)enddodo k=1,kepv(1:ie,1:je,k)=1.E6*9.80665*(> -(-dudy(1:ie,1:je,k)+dvdx(1:ie,1:je,k)> +f(1:ie,1:je))*dthdp(1:ie,1:je,k)> -(dudp(1:ie,1:je,k)*dthdy(1:ie,1:je,k)> -dvdp(1:ie,1:je,k)*dthdx(1:ie,1:je,k)))enddoc interpolate poles from neighbouring points on every levelif (vmax(2).gt.89.5) thendo k=1,kepvpole=0.do i=1,iepvpole=pvpole+pv(i,je-1,k)enddopvpole=pvpole/real(ie)do i=1,iepv(i,je,k)=pvpoleenddoenddoendifif (vmin(2).lt.-89.5) thendo k=1,kepvpole=0.do i=1,iepvpole=pvpole+pv(i,2,k)enddopvpole=pvpole/real(ie)do i=1,iepv(i,1,k)=pvpoleenddoenddoendifIF (ALLOCATED(dspdx)) DEALLOCATE(dspdx)IF (ALLOCATED(dspdy)) DEALLOCATE(dspdy)IF (ALLOCATED(dthdp)) DEALLOCATE(dthdp)IF (ALLOCATED(dudp)) DEALLOCATE(dudp)IF (ALLOCATED(dvdp)) DEALLOCATE(dvdp)IF (ALLOCATED(dvdx)) DEALLOCATE(dvdx)IF (ALLOCATED(dudy)) DEALLOCATE(dudy)IF (ALLOCATED(dthdx)) DEALLOCATE(dthdx)IF (ALLOCATED(dthdy)) DEALLOCATE(dthdy)IF (ALLOCATED(uu2)) DEALLOCATE(uu2)endsubroutine potvort_i(pv,uu,vv,p,cl,f,ie,je,ke,ak,> vmin,vmax,mdv)C =================================================C calculate isentropic PVC !! calculating PV makes sense onlyC when theta levels are very close to each others !!C with misdat (mdv) treatmentc argument declarationinteger ie,je,kereal pv(ie,je,ke),uu(ie,je,ke),vv(ie,je,ke),> p(ie,je,ke),cl(ie,je),f(ie,je)real ak(ke),vmin(4),vmax(4),mdvreal pvpolec variable declarationREAL,ALLOCATABLE, DIMENSION (:,:,:) :: dpdthREAL,ALLOCATABLE, DIMENSION (:,:,:) :: uu2,dvdx,dudyinteger i,j,k,statallocate(dpdth(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort_i: error allocating dpdth'allocate(dvdx(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort_i: error allocating dvdx'allocate(dudy(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort_i: error allocating dudy'allocate(uu2(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort_i: error allocating uu2'c calculate dp/dthc k=1do i=1,iedo j=1,jeif ((p(i,j,2).eq.mdv).or.(p(i,j,1).eq.mdv)) thendpdth(i,j,1)=mdvelsedpdth(i,j,1)=100.*(p(i,j,2)-p(i,j,1))/(ak(2)-ak(1))endifenddoenddoc k=2,ke-1do i=1,iedo j=1,jedo k=2,ke-1if ((p(i,j,k+1).eq.mdv).or.(p(i,j,k-1).eq.mdv)) thendpdth(i,j,k)=mdvelsedpdth(i,j,k)=100.*(p(i,j,k+1)-p(i,j,k-1))/(ak(k+1)-ak(k-1))endifenddoenddoenddoc k=kedo i=1,iedo j=1,jeif ((p(i,j,ke).eq.mdv).or.(p(i,j,ke-1).eq.mdv)) thendpdth(i,j,ke)=mdvelsedpdth(i,j,ke)=100.*(p(i,j,ke)-p(i,j,ke-1))/(ak(ke)-ak(ke-1))endifenddoenddocall ddh2m(vv,dvdx,cl,'X',ie,je,ke,vmin,vmax,mdv)c conversion of uu for spheric coordinates (uu*cos(phi))do i=1,iedo j=1,jedo k=1,keif (uu(i,j,k).ne.mdv) thenuu2(i,j,k)=uu(i,j,k)*cl(i,j)elseuu2(i,j,k)=mdvendifenddoenddoenddocall ddh2m(uu2,dudy,cl,'Y',ie,je,ke,vmin,vmax,mdv)c conversion of dudy for spheric coordinates (dudy/cos(phi))do i=1,iedo j=1,jedo k=1,keif (dudy(i,j,k).ne.mdv) thendudy(i,j,k)=dudy(i,j,k)/cl(i,j)elsedudy(i,j,k)=mdvendifenddoenddoenddodo i=1,iedo j=1,jedo k=1,keif ((dudy(i,j,k).eq.mdv).or.(dvdx(i,j,k).eq.mdv).or.> (dpdth(i,j,k).eq.mdv).or.(dpdth(i,j,k).eq.0.)) thenpv(i,j,k)=mdvelsepv(i,j,k)=-1.E6*9.80665*> (-dudy(i,j,k)+dvdx(i,j,k)+f(i,j))/dpdth(i,j,k)endifenddoenddoenddoc interpolate poles from neighbouring points on every levelif (vmax(2).eq.90.) thendo k=1,kepvpole=0.do i=1,iepvpole=pvpole+pv(i,je-1,k)enddopvpole=pvpole/real(ie)do i=1,iepv(i,je,k)=pvpoleenddoenddoendifif (vmin(2).eq.-90.) thendo k=1,kepvpole=0.do i=1,iepvpole=pvpole+pv(i,2,k)enddopvpole=pvpole/real(ie)do i=1,iepv(i,1,k)=pvpoleenddoenddoendifIF (ALLOCATED(dpdth)) DEALLOCATE(dpdth)IF (ALLOCATED(dvdx)) DEALLOCATE(dvdx)IF (ALLOCATED(dudy)) DEALLOCATE(dudy)IF (ALLOCATED(uu2)) DEALLOCATE(uu2)endsubroutine gpotvort_i2(dpv,uu,vv,p,cl,f,ie,je,ke,ak,> vmin,vmax,mdv)C =================================================C calculate isentropic PV gradient lengthC GPV=SQRT(D[PV:X]^2+D[PV:Y]^2)C !! it calculates PV with potvort_i, what makes sense onlyC when theta levels are very close to each others !!C with misdat (mdv) treatmentc argument declarationinteger ie,je,kereal dpv(ie,je,ke),uu(ie,je,ke),vv(ie,je,ke),> p(ie,je,ke),cl(ie,je),f(ie,je)real ak(ke),vmin(4),vmax(4),mdvc variable declarationREAL,ALLOCATABLE, DIMENSION (:,:,:) :: pvREAL,ALLOCATABLE, DIMENSION (:,:) :: dpvdx,dpvdyinteger statreal dpvpoleallocate(pv(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'gpotvort_i2: error allocating pv'allocate(dpvdx(ie,je),STAT=stat)if (stat.ne.0) print*,'gpotvort_i2: error allocating dpvdx'allocate(dpvdy(ie,je),STAT=stat)if (stat.ne.0) print*,'gpotvort_i2: error allocating dpvdy'call potvort_i(pv,uu,vv,p,cl,f,ie,je,ke,ak,> vmin,vmax,mdv)do k=1,kecall ddh2m(pv(1,1,k),dpvdx,cl,'X',ie,je,1,vmin,vmax,mdv)call ddh2m(pv(1,1,k),dpvdy,cl,'Y',ie,je,1,vmin,vmax,mdv)do j=1,jedo i=1,ieif ((dpvdx(i,j).ne.mdv).and.(dpvdy(i,j).ne.mdv)) thendpv(i,j,k)=1.E6*sqrt(dpvdx(i,j)**2.+dpvdy(i,j)**2.)elsedpv(i,j,k)=mdvendifenddoenddoenddoIF (ALLOCATED(pv)) DEALLOCATE(pv)IF (ALLOCATED(dpvdx)) DEALLOCATE(dpvdx)IF (ALLOCATED(dpvdy)) DEALLOCATE(dpvdy)endsubroutine gpotvort_i(dpv,pv,cl,ie,je,ke,> vmin,vmax,mdv)C =================================================C calculate isentropic PV gradient lengthC GPV=SQRT(D[PV:X]^2+D[PV:Y]^2)C !! it uses PV as input, best calculated withC potvort (with data from P file) !!C with misdat (mdv) treatmentC 990201 mark linigerc argument declarationinteger ie,je,kereal dpv(ie,je,ke),pv(ie,je,ke),cl(ie,je)real vmin(4),vmax(4),mdvc variable declarationREAL,ALLOCATABLE, DIMENSION (:,:) :: dpvdx,dpvdyinteger statallocate(dpvdx(ie,je),STAT=stat)if (stat.ne.0) print*,'gpotvort_i: error allocating dpvdx'allocate(dpvdy(ie,je),STAT=stat)if (stat.ne.0) print*,'gpotvort_i: error allocating dpvdy'do k=1,kecall ddh2m(pv(1,1,k),dpvdx,cl,'X',ie,je,1,vmin,vmax,mdv)call ddh2m(pv(1,1,k),dpvdy,cl,'Y',ie,je,1,vmin,vmax,mdv)do j=1,jedo i=1,ieif ((dpvdx(i,j).ne.mdv).and.(dpvdy(i,j).ne.mdv)) thendpv(i,j,k)=1.E6*sqrt(dpvdx(i,j)**2.+dpvdy(i,j)**2.)elsedpv(i,j,k)=mdvendifenddoenddoenddoIF (ALLOCATED(dpvdx)) DEALLOCATE(dpvdx)IF (ALLOCATED(dpvdy)) DEALLOCATE(dpvdy)endsubroutine divqu(var,uu,vv,qq,cl,ie,je,ke,vmin,vmax,mdv)C ========================================================C calculate divergence DIV=D[U:X]+D[V*cos(phi):Y]/cos(phi)C with misdat (mdv) and pole treatmentc argument declarationinteger ie,je,kereal,intent(OUT) :: var(ie,je,ke)real,intent(IN) :: uu(ie,je,ke),vv(ie,je,ke),qq(ie,je,ke),> cl(ie,je)real vmin(4),vmax(4),mdvc variable declarationREAL,ALLOCATABLE, DIMENSION (:,:) :: qv2,dqudx,dqvdyREAL,ALLOCATABLE, DIMENSION (:,:,:) :: qu,qvinteger statreal poleallocate(qv2(ie,je),STAT=stat)if (stat.ne.0) print*,'divqu: error allocating qv2'allocate(dqudx(ie,je),STAT=stat)if (stat.ne.0) print*,'divqu: error allocating dqudx'allocate(dqvdy(ie,je),STAT=stat)if (stat.ne.0) print*,'divqu: error allocating dqvdy'allocate(qu(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'divqu: error allocating qu'allocate(qv(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'divqu: error allocating qv'do k=1,kedo j=1,jedo i=1,iequ(i,j,k)=qq(i,j,k)*uu(i,j,k)qv(i,j,k)=qq(i,j,k)*vv(i,j,k)enddoenddoenddodo k=1,kecall ddh2m(qu(1,1,k),dqudx,cl,'X',ie,je,1,vmin,vmax,mdv)c conversion of qv for spheric coordinates (qv*cos(phi))do j=1,jedo i=1,ieif (qv(i,j,k).ne.mdv) thenqv2(i,j)=qv(i,j,k)*cl(i,j)elseqv2(i,k)=mdvendifenddoenddocall ddh2m(qv2,dqvdy,cl,'Y',ie,je,1,vmin,vmax,mdv)c conversion of dqvdy for spheric coordinates (dqvdy/cos(phi))do j=1,jedo i=1,ieif ((dqudx(i,j).ne.mdv).and.(dqvdy(i,j).ne.mdv)) thenvar(i,j,k)=1.E6*(dqudx(i,j)+dqvdy(i,j)/cl(i,j))elsevar(i,j,k)=mdvendifenddoenddoenddoc interpolate poles from neighbouring points on every levelif (vmax(2).eq.90.) thendo k=1,kepole=0.counter=0.do i=1,ieif (var(i,je-1,k).ne.mdv) thencounter=counter+1pole=pole+var(i,je-1,k)endifenddoif (counter.ne.0) thenpole=pole/counterelsepole=mdvendifdo i=1,ievar(i,je,k)=poleenddoenddoendifif (vmin(2).eq.-90.) thendo k=1,kepole=0.counter=0.do i=1,ieif (var(i,2,k).ne.mdv) thencounter=counter+1pole=pole+var(i,2,k)endifenddoif (counter.ne.0) thenpole=pole/counterelsepole=mdvendifdo i=1,ievar(i,1,k)=poleenddoenddoendifIF (ALLOCATED(dqudx)) DEALLOCATE(dqudx)IF (ALLOCATED(dqvdy)) DEALLOCATE(dqvdy)IF (ALLOCATED(qu)) DEALLOCATE(qu)IF (ALLOCATED(qv)) DEALLOCATE(qv)IF (ALLOCATED(qv2)) DEALLOCATE(qv2)endsubroutine div_i(var,uu,vv,cl,ie,je,ke,vmin,vmax,mdv)C =====================================================C calculate divergence DIV=D[U:X]+D[V*cos(phi):Y]/cos(phi)C with misdat (mdv) and pole treatmentC 990201 mark linigerc argument declarationinteger ie,je,kereal,intent(OUT) :: var(ie,je,ke)real,intent(IN) :: uu(ie,je,ke),vv(ie,je,ke),cl(ie,je)real vmin(4),vmax(4),mdvc variable declarationREAL,ALLOCATABLE, DIMENSION (:,:) :: vv2,dudx,dvdyinteger statreal poleallocate(vv2(ie,je),STAT=stat)if (stat.ne.0) print*,'div_i: error allocating vv2'allocate(dudx(ie,je),STAT=stat)if (stat.ne.0) print*,'div_i: error allocating dudx'allocate(dvdy(ie,je),STAT=stat)if (stat.ne.0) print*,'div_i: error allocating dvdy'do k=1,kecall ddh2m(uu(1,1,k),dudx,cl,'X',ie,je,1,vmin,vmax,mdv)c conversion of vv for spheric coordinates (vv*cos(phi))do j=1,jedo i=1,ieif (vv(i,j,k).ne.mdv) thenvv2(i,j)=vv(i,j,k)*cl(i,j)elsevv2(i,k)=mdvendifenddoenddocall ddh2m(vv2,dvdy,cl,'Y',ie,je,1,vmin,vmax,mdv)c conversion of dvdy for spheric coordinates (dvdy/cos(phi))do j=1,jedo i=1,ieif ((dudx(i,j).ne.mdv).and.(dvdy(i,j).ne.mdv)) thenvar(i,j,k)=1.E6*(dudx(i,j)+dvdy(i,j)/cl(i,j))elsevar(i,j,k)=mdvendifenddoenddoenddoc interpolate poles from neighbouring points on every levelif (vmax(2).eq.90.) thendo k=1,kepole=0.counter=0.do i=1,ieif (var(i,je-1,k).ne.mdv) thencounter=counter+1pole=pole+var(i,je-1,k)endifenddoif (counter.ne.0) thenpole=pole/counterelsepole=mdvendifdo i=1,ievar(i,je,k)=poleenddoenddoendifif (vmin(2).eq.-90.) thendo k=1,kepole=0.counter=0.do i=1,ieif (var(i,2,k).ne.mdv) thencounter=counter+1pole=pole+var(i,2,k)endifenddoif (counter.ne.0) thenpole=pole/counterelsepole=mdvendifdo i=1,ievar(i,1,k)=poleenddoenddoendifIF (ALLOCATED(dudx)) DEALLOCATE(dudx)IF (ALLOCATED(vv2)) DEALLOCATE(vv2)IF (ALLOCATED(dvdy)) DEALLOCATE(dvdy)endsubroutine def_i(var,uu,vv,cl,ie,je,ke,vmin,vmax,mdv)C =====================================================C calculate deformation (strain) DEF= 1e6*SQRT(C ( D[V:X]+D[U*cos(phi):Y]/cos(phi) )^2C ( D[U:X]-D[V*cos(phi):Y]/cos(phi) )^2 )C with misdat (mdv) treatmentC 990501 mark linigerc argument declarationinteger ie,je,kereal,intent(OUT) :: var(ie,je,ke)real,intent(IN) :: uu(ie,je,ke),vv(ie,je,ke),cl(ie,je)real vmin(4),vmax(4),mdvc variable declarationREAL,ALLOCATABLE, DIMENSION (:,:) :: uu2,vv2,dudx,dvdx,dudy,dvdyinteger statreal poleallocate(uu2(ie,je),STAT=stat)if (stat.ne.0) print*,'dev_i: error allocating uu2'allocate(vv2(ie,je),STAT=stat)if (stat.ne.0) print*,'dev_i: error allocating vv2'allocate(dudx(ie,je),STAT=stat)if (stat.ne.0) print*,'dev_i: error allocating dudx'allocate(dvdx(ie,je),STAT=stat)if (stat.ne.0) print*,'dev_i: error allocating dvdx'allocate(dudy(ie,je),STAT=stat)if (stat.ne.0) print*,'dev_i: error allocating dudy'allocate(dvdy(ie,je),STAT=stat)if (stat.ne.0) print*,'dev_i: error allocating dvdy'do k=1,kecall ddh2m(uu(1,1,k),dudx,cl,'X',ie,je,1,vmin,vmax,mdv)call ddh2m(vv(1,1,k),dvdx,cl,'X',ie,je,1,vmin,vmax,mdv)c conversion of uu/vv for spheric coordinates (??*cos(phi))do j=1,jedo i=1,ieif (uu(i,j,k).ne.mdv) thenuu2(i,j)=uu(i,j,k)*cl(i,j)elseuu2(i,j)=mdvendifif (vv(i,j,k).ne.mdv) thenvv2(i,j)=vv(i,j,k)*cl(i,j)elsevv2(i,j)=mdvendifenddoenddocall ddh2m(uu2,dudy,cl,'Y',ie,je,1,vmin,vmax,mdv)call ddh2m(vv2,dvdy,cl,'Y',ie,je,1,vmin,vmax,mdv)c conversion of d?dy for spheric coordinates (d?dy/cos(phi))do j=1,jedo i=1,ieif ((dudx(i,j).ne.mdv).and.> (dudy(i,j).ne.mdv).and.> (dvdx(i,j).ne.mdv).and.> (dvdy(i,j).ne.mdv)) thenvar(i,j,k)=1.e6*sqrt(> (dvdx(i,j)+dudy(i,j)/cl(i,j))> *(dvdx(i,j)+dudy(i,j)/cl(i,j))> +(dudx(i,j)-dvdy(i,j)/cl(i,j))> *(dudx(i,j)-dvdy(i,j)/cl(i,j)))elsevar(i,j,k)=mdvendifenddoenddoenddoc interpolate poles from neighbouring points on every levelif (vmax(2).eq.90.) thendo k=1,kepole=0.do i=1,iepole=pole+var(i,je-1,k)enddopole=pole/real(ie)do i=1,ievar(i,je,k)=poleenddoenddoendifif (vmin(2).eq.-90.) thendo k=1,kepole=0.do i=1,iepole=pole+var(i,2,k)enddopole=pole/real(ie)do i=1,ievar(i,1,k)=poleenddoenddoendifIF (ALLOCATED(uu2)) DEALLOCATE(uu2)IF (ALLOCATED(vv2)) DEALLOCATE(vv2)IF (ALLOCATED(dudx)) DEALLOCATE(dudx)IF (ALLOCATED(dudy)) DEALLOCATE(dudy)IF (ALLOCATED(dvdx)) DEALLOCATE(dvdx)IF (ALLOCATED(dvdy)) DEALLOCATE(dvdy)endsubroutine defn_i(var,uu,vv,cl,ie,je,ke,vmin,vmax,mdv)C =====================================================C calculate normal deformation (normal strain)C DEFN= 1e6*( D[U:X]-D[V*cos(phi):Y]/cos(phi) )C with misdat (mdv) treatmentC 990727 mark linigerc argument declarationinteger ie,je,kereal,intent(OUT) :: var(ie,je,ke)real,intent(IN) :: uu(ie,je,ke),vv(ie,je,ke),cl(ie,je)real vmin(4),vmax(4),mdvc variable declarationREAL,ALLOCATABLE, DIMENSION (:,:) :: vv2,dudx,dvdyinteger statreal poleallocate(vv2(ie,je),STAT=stat)if (stat.ne.0) print*,'dev_i: error allocating vv2'allocate(dudx(ie,je),STAT=stat)if (stat.ne.0) print*,'dev_i: error allocating dudx'allocate(dvdy(ie,je),STAT=stat)if (stat.ne.0) print*,'dev_i: error allocating dvdy'do k=1,kecall ddh2m(uu(1,1,k),dudx,cl,'X',ie,je,1,vmin,vmax,mdv)c conversion of uu/vv for spheric coordinates (??*cos(phi))do j=1,jedo i=1,ieif (vv(i,j,k).ne.mdv) thenvv2(i,j)=vv(i,j,k)*cl(i,j)elsevv2(i,j)=mdvendifenddoenddocall ddh2m(vv2,dvdy,cl,'Y',ie,je,1,vmin,vmax,mdv)c final calculationdo j=1,jedo i=1,ieif ((dudx(i,j).ne.mdv).and.> (dvdy(i,j).ne.mdv)) thenvar(i,j,k)=1.e6*(dudx(i,j)-dvdy(i,j)/cl(i,j))elsevar(i,j,k)=mdvendifenddoenddoenddoc interpolate poles from neighbouring points on every levelif (vmax(2).eq.90.) thendo k=1,kepole=0.do i=1,iepole=pole+var(i,je-1,k)enddopole=pole/real(ie)do i=1,ievar(i,je,k)=poleenddoenddoendifif (vmin(2).eq.-90.) thendo k=1,kepole=0.do i=1,iepole=pole+var(i,2,k)enddopole=pole/real(ie)do i=1,ievar(i,1,k)=poleenddoenddoendifIF (ALLOCATED(vv2)) DEALLOCATE(vv2)IF (ALLOCATED(dudx)) DEALLOCATE(dudx)IF (ALLOCATED(dvdy)) DEALLOCATE(dvdy)endsubroutine defs_i(var,uu,vv,cl,ie,je,ke,vmin,vmax,mdv)C =====================================================C calculate shear deformation (shear strain)C DEFS= 1e6*( D[V:X]+D[U*cos(phi):Y]/cos(phi))C with misdat (mdv) treatmentC 990727 mark linigerc argument declarationinteger ie,je,kereal,intent(OUT) :: var(ie,je,ke)real,intent(IN) :: uu(ie,je,ke),vv(ie,je,ke),cl(ie,je)real vmin(4),vmax(4),mdvc variable declarationREAL,ALLOCATABLE, DIMENSION (:,:) :: uu2,dvdx,dudyinteger statreal poleallocate(uu2(ie,je),STAT=stat)if (stat.ne.0) print*,'dev_i: error allocating uu2'allocate(dvdx(ie,je),STAT=stat)if (stat.ne.0) print*,'dev_i: error allocating dvdx'allocate(dudy(ie,je),STAT=stat)if (stat.ne.0) print*,'dev_i: error allocating dudy'do k=1,kecall ddh2m(vv(1,1,k),dvdx,cl,'X',ie,je,1,vmin,vmax,mdv)c conversion of uu/vv for spheric coordinates (??*cos(phi))do j=1,jedo i=1,ieif (uu(i,j,k).ne.mdv) thenuu2(i,j)=uu(i,j,k)*cl(i,j)elseuu2(i,j)=mdvendifenddoenddocall ddh2m(uu2,dudy,cl,'Y',ie,je,1,vmin,vmax,mdv)c conversion of d?dy for spheric coordinates (d?dy/cos(phi))do j=1,jedo i=1,ieif ((dudy(i,j).ne.mdv).and.> (dvdx(i,j).ne.mdv)) thenvar(i,j,k)=1.e6*(dvdx(i,j)+dudy(i,j)/cl(i,j))elsevar(i,j,k)=mdvendifenddoenddoenddoc interpolate poles from neighbouring points on every levelif (vmax(2).eq.90.) thendo k=1,kepole=0.do i=1,iepole=pole+var(i,je-1,k)enddopole=pole/real(ie)do i=1,ievar(i,je,k)=poleenddoenddoendifif (vmin(2).eq.-90.) thendo k=1,kepole=0.do i=1,iepole=pole+var(i,2,k)enddopole=pole/real(ie)do i=1,ievar(i,1,k)=poleenddoenddoendifIF (ALLOCATED(uu2)) DEALLOCATE(uu2)IF (ALLOCATED(dudy)) DEALLOCATE(dudy)IF (ALLOCATED(dvdx)) DEALLOCATE(dvdx)endsubroutine okuboweiss_i(var,uu,vv,cl,ie,je,ke,vmin,vmax,mdv)C =====================================================C calculate okubo-weiss parameter= DEF*DEF-VORT*VORTC with misdat (mdv) treatmentC 990720 mark linigerc argument declarationinteger ie,je,kereal,intent(OUT) :: var(ie,je,ke)real,intent(IN) :: uu(ie,je,ke),vv(ie,je,ke),cl(ie,je)real vmin(4),vmax(4),mdvc variable declarationREAL,ALLOCATABLE, DIMENSION (:,:,:) :: def,vortinteger statreal poleallocate(def(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'okuboweiss_i: error allocating def'allocate(vort(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'okuboweiss_i: error allocating vort'call def_i(def,uu,vv,cl,ie,je,ke,vmin,vmax,mdv)call vort_i(vort,uu,vv,cl,ie,je,ke,vmin,vmax,mdv)do k=1,kedo j=1,jedo i=1,ieif ((def(i,j,k).ne.mdv).and.> (vort(i,j,k).ne.mdv)) thenvar(i,j,k)=def(i,j,k)*def(i,j,k)> -vort(i,j,k)*vort(i,j,k)*1e4elsevar(i,j,k)=mdvendifenddoenddoenddoIF (ALLOCATED(def)) DEALLOCATE(def)IF (ALLOCATED(vort)) DEALLOCATE(vort)endsubroutine rich(ri,sp,uu,vv,th,ie,je,ke,ak,bk)C ==============================================C calculate Richardson numberreal R,p0,kappadata R,p0,kappa /287.,100000.,0.286/c argument declarationinteger ie,je,kereal ri(ie,je,ke),uu(ie,je,ke),vv(ie,je,ke),> th(ie,je,ke),sp(ie,je)real ak(ke),bk(ke)c variable declarationREAL,ALLOCATABLE, DIMENSION (:,:,:) :: dthdp,dudp,dvdpinteger i,j,k,statallocate(dthdp(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort: error allocating dthdp'* allocate(vel(ie,je,ke),STAT=stat)* if (stat.ne.0) print*,'potvort: error allocating vel'allocate(dudp(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort: error allocating dudp'allocate(dvdp(ie,je,ke),STAT=stat)if (stat.ne.0) print*,'potvort: error allocating dvdp'* vel=sqrt(uu**2+vv**2)call ddp(th,dthdp,sp,ie,je,ke,ak,bk)* call ddp(vel,dveldp,sp,ie,je,ke,ak,bk)call ddp(uu,dudp,sp,ie,je,ke,ak,bk)call ddp(vv,dvdp,sp,ie,je,ke,ak,bk)do k=1,kedo j=1,jedo i=1,iepval=(ak(k)+bk(k)*sp(i,j))*100.ri(i,j,k)=-R*(p0**(-kappa))*(pval**(kappa-1.))*> dthdp(i,j,k)/(dudp(i,j,k)**2.+dvdp(i,j,k)**2.)* ri(i,j,k)=-R*(p0**(-kappa))*(pval**(kappa-1.))** > dthdp(i,j,k)/(dveldp(i,j,k)**2.)enddoenddoenddoIF (ALLOCATED(dthdp)) DEALLOCATE(dthdp)IF (ALLOCATED(dudp)) DEALLOCATE(dudp)IF (ALLOCATED(dvdp)) DEALLOCATE(dvdp)endsubroutine ddp(a,d,sp,ie,je,ke,ak,bk)c-----------------------------------------------------------------------c Purpose: VERTICAL DERIVATIVEc Compute the vertical derivative without missing data checking.c The derivative is taken from array 'a' in the direction of 'P'.c The result is stored in array 'd'.c 3 point weighted centered differencing is used.c The vertical level-structure of the data is of the formc p=ak(k)+bk(k)*ps.c-----------------------------------------------------------------------c declaration of argumentsinteger ie,je,kereal a(ie,je,ke),d(ie,je,ke),sp(ie,je)c variable declarationinteger i,j,kreal dpu,dpl,quot,fac,psrfreal ak(ke),bk(ke)c specify scaling factor associated with derivation with respectc to pressurefac=0.01c compute vertical 3 point derivativec ---------------------------c 3-point vertical derivativec ---------------------------do j=1,jedo i=1,iec get surface pressure at current grid-pointpsrf=sp(i,j)c points at k=1dpu=(ak(1)+bk(1)*psrf)-(ak(2)+bk(2)*psrf)d(i,j,1)=(a(i,j,1)-a(i,j,2))*fac/dpuc points at 1<k<kedo k=2,ke-1dpu=(ak(k)+bk(k)*psrf)-(ak(k+1)+bk(k+1)*psrf)dpl=(ak(k-1)+bk(k-1)*psrf)-(ak(k)+bk(k)*psrf)quot=dpu/dpld(i,j,k)=(quot*(a(i,j,k-1)-a(i,j,k))& +1./quot*(a(i,j,k)-a(i,j,k+1)))*fac/(dpu+dpl)enddoc points at k=kedpl=(ak(ke-1)+bk(ke-1)*psrf)-(ak(ke)+bk(ke)*psrf)d(i,j,ke)=(a(i,j,ke-1)-a(i,j,ke))*fac/dplenddoenddoendsubroutine ddh3(a,d,ps,dps,cl,dir,ie,je,ke,datmin,datmax,ak,bk)c-----------------------------------------------------------------------c Purpose: HORIZONTAL DERIVATIVE ON PRESSURE-SURFACES WITHOUTc MISSING DATAc The derivative is taken from array 'a' in the direction of 'dir',c where 'dir' is either 'X','Y'. The result is stored in array 'd'.c The routine accounts for derivatives at the pole and periodicc boundaries in the longitudinal direction (depending onc the value of datmin, datmax). If the data-set does not reach toc the pole, a one-sided derivative is taken. Pole-treatment is onlyc carried out if the data-set covers 360 deg in longitude, and itc requires that ie=4*ii+1, where ii is an integer.c History:c Daniel Luethic-----------------------------------------------------------------------c declaration of argumentsinteger ie,je,kereal a(ie,je,ke),d(ie,je,ke),cl(ie,je)real ps(ie,je),dps(ie,je)real datmin(4),datmax(4)character*(*) dirc variable declarationinteger i,j,kreal ak(ke),bk(ke),as(500),bs(500)c compute vertical derivatives of ak's and bk'sdo k=2,ke-1as(k)=(ak(k-1)-ak(k+1))/2.bs(k)=(bk(k-1)-bk(k+1))/2.enddoas(1 )=ak(1)-ak(2)bs(1 )=bk(1)-bk(2)as(ke)=ak(ke-1)-ak(ke)bs(ke)=bk(ke-1)-bk(ke)c compute horizontal derivatives on sigma surfacescall ddh2(a,d,cl,dir,ie,je,ke,datmin,datmax)c apply correction for horizontal derivative on p-surfacesdo j=1,jedo i=1,iedo k=2,ke-1d(i,j,k)=d(i,j,k)+bk(k)*dps(i,j)/2./(as(k)+& bs(k)*ps(i,j))*(a(i,j,k+1)-a(i,j,k-1))enddok=1d(i,j,k)=d(i,j,k)+bk(k)*dps(i,j)/(as(k)+& bs(k)*ps(i,j))*(a(i,j,k+1)-a(i,j,k))k=ked(i,j,k)=d(i,j,k)+bk(k)*dps(i,j)/(as(k)+& bs(k)*ps(i,j))*(a(i,j,k)-a(i,j,k-1))enddoenddoendsubroutine ddh2(a,d,cl,dir,ie,je,ke,datmin,datmax)c-----------------------------------------------------------------------c Purpose: HORIZONTAL DERIVATIVE ON DATA-SURFACES WITHOUT MISSINGC DATAc Compute the horizontal derivative without missing data checking.c The derivative is taken from array 'a' in the direction of 'dir',c where 'dir' is either 'X','Y'. The result is stored in array 'd'.c The routine accounts for derivatives at the pole and periodicc boundaries in the longitudinal direction (depending onc the value of datmin, datmax). If the data-set does not reach toc the pole, a one-sided derivative is taken. Pole-treatment is onlyc carried out if the data-set covers 360 deg in longitude, and itc requires that ie=4*ii+1, where ii is an integer.c-----------------------------------------------------------------------c declaration of argumentsinteger ie,je,kereal a(ie,je,ke),d(ie,je,ke),cl(ie,je)real datmin(4),datmax(4)character*(*) dirc local variable declarationinteger i,j,k,ip1,im1,jp1,jm1,ip,im,j1,j2real dlat,dlon,coslat,dx,dy,dxr,dyrinteger northpl,southpl,lonperc rerd and circ are the mean radius and diameter of the earth inc meterreal rerd,circ,pidata rerd,circ,pi /6.37e6,4.e7,3.141592654/c compute flags for pole and periodic treatmentsouthpl=0northpl=0lonper =0j1=1j2=jeif (abs(datmax(1)-datmin(1)-360.).lt.1.e-3) thenlonper=1if (abs(datmin(2)+90.).lt.1.e-3) thensouthpl=1j1=2endifif (abs(datmax(2)-90.).lt.1.e-3) thennorthpl=1j2=je-1endifendifdlon=((datmax(1)-datmin(1))/float(ie-1)) *pi/180.dlat=((datmax(2)-datmin(2))/float(je-1)) *pi/180.c print *,'Computing derivative ',dir(1:1),c & ' of an array dimensioned ',ie,je,keif (dir(1:1).eq.'X') thenc -----------------------------c derivation in the x-directionc -----------------------------do k=1,kec do gridpoints at j1<=j<=j2do j=j1,j2coslat=cl(1,j)c do regular gridpoints at 1<i<ie, 1<j<jedx =rerd*coslat*dlondxr=1./(2.*dx)do i=2,ie-1ip1=i+1im1=i-1d(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))enddo ! i-loopc completed regular gridpoints at 1<i<ie, 1<j<jec do gridpoints at i=1, i=ie, 1<j<jeif (lonper.eq.1) thenc use periodic boundariesi=1ip1=2im1=ie-1d(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))d(ie,j,k)=d(1,j,k)elsec use one-sided derivativesdxr=1./dxi=1ip1=2im1=1d(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))i=ieip1=ieim1=ie-1d(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))endifc completed gridpoints at i=1, i=ie, j1<=j<=j2enddo ! j-loopc completed gridpoints at 1<j<jec do gridpoints at j=jeif (northpl.eq.1) thenc for these gridpoints, the derivative in the x-direction is ac derivative in the y-direction at another pole-gridpointdy =rerd*dlatdyr=1./(2.*dy)j=jejp1=je-1jm1=je-1do i=1,ieip=mod(i-1+ (ie-1)/4,ie)+1im=mod(i-1+3*(ie-1)/4,ie)+1d(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))enddo ! i-loopc completed gridpoints at j=jeendifc do gridpoints at j=1if (southpl.eq.1) thendy =rerd*dlatdyr=1./(2.*dy)j=1jp1=2jm1=2do i=1,ieip=mod(i-1+ (ie-1)/4,ie)+1im=mod(i-1+3*(ie-1)/4,ie)+1d(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))enddo ! i-loopendifc completed gridpoints at j=1enddo ! k-loopelse if (dir(1:1).eq.'Y') thenc -----------------------------c derivation in the y-directionc -----------------------------dy =dlat*rerddyr=1./(2.*dy)do k=1,kedo i=1,iec do regular gridpointsdo j=2,je-1jp1=j+1jm1=j-1d(i,j,k)=dyr*(a(i,jp1,k)-a(i,jm1,k))enddoc do gridpoints at j=jeif (northpl.eq.1) thenc pole-treatmentj=jejm1=j-1jp1=j-1ip=mod(i-1+(ie-1)/2,ie)+1im=id(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))elsec one-sided derivativej=jejm1=j-1jp1=jd(i,j,k)=2.*dyr*(a(i,jp1,k)-a(i,jm1,k))endifc completed gridpoints at j=jec do gridpoints at j=1if (southpl.eq.1) thenc pole-treatmentj=1jm1=2jp1=2ip=iim=mod(i-1+(ie-1)/2,ie)+1d(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))elsec one-sided derivativej=1jm1=1jp1=2d(i,j,k)=2.*dyr*(a(i,jp1,k)-a(i,jm1,k))endifc completed gridpoints at j=1enddoenddoendifendsubroutine ddh2m(a,d,cl,dir,ie,je,ke,datmin,datmax,mdv)c-----------------------------------------------------------------------c Purpose: HORIZONTAL DERIVATIVE ON DATA-SURFACES WITH MISSINGC DATAc Compute the horizontal derivative with missing data checking.c The derivative is taken from array 'a' in the direction of 'dir',c where 'dir' is either 'X','Y'. The result is stored in array 'd'.c The routine accounts for derivatives at the pole and periodicc boundaries in the longitudinal direction (depending onc the value of datmin, datmax). If the data-set does not reach toc the pole, a one-sided derivative is taken. Pole-treatment is onlyc carried out if the data-set covers 360 deg in longitude, and itc requires that ie=4*ii+1, where ii is an integer.c mdv is the misdat parameter (real)c-----------------------------------------------------------------------c declaration of argumentsinteger ie,je,kereal a(ie,je,ke),d(ie,je,ke),cl(ie,je)real datmin(4),datmax(4),mdvcharacter*(*) dirc local variable declarationinteger i,j,k,ip1,im1,jp1,jm1,ip,im,j1,j2real dlat,dlon,coslat,dx,dy,dxr,dyrinteger northpl,southpl,lonperc rerd and circ are the mean radius and diameter of the earth inc meterreal rerd,circ,pidata rerd,circ,pi /6.37e6,4.e7,3.141592654/c compute flags for pole and periodic treatmentsouthpl=0northpl=0lonper =0j1=1j2=jeif (abs(datmax(1)-datmin(1)-360.).lt.1.e-3) thenlonper=1if (abs(datmin(2)+90.).lt.1.e-3) thensouthpl=1j1=2endifif (abs(datmax(2)-90.).lt.1.e-3) thennorthpl=1j2=je-1endifendifdlon=((datmax(1)-datmin(1))/float(ie-1)) *pi/180.dlat=((datmax(2)-datmin(2))/float(je-1)) *pi/180.c print *,'Computing derivative ',dir(1:1),c & ' of an array dimensioned ',ie,je,keif (dir(1:1).eq.'X') thenc -----------------------------c derivation in the x-directionc -----------------------------do k=1,kec do gridpoints at j1<=j<=j2do j=j1,j2coslat=cl(1,j)c do regular gridpoints at 1<i<ie, 1<j<jedx =rerd*coslat*dlondxr=1./(2.*dx)do i=2,ie-1ip1=i+1im1=i-1if ((a(ip1,j,k).ne.mdv).and.(a(im1,j,k).ne.mdv)) thend(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))elsed(i,j,k)=mdvendifenddo ! i-loopc completed regular gridpoints at 1<i<ie, 1<j<jec do gridpoints at i=1, i=ie, 1<j<jeif (lonper.eq.1) thenc use periodic boundariesi=1ip1=2im1=ie-1if ((a(ip1,j,k).ne.mdv).and.(a(im1,j,k).ne.mdv)) thend(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))elsed(i,j,k)=mdvendifd(ie,j,k)=d(1,j,k)elsec use one-sided derivativesdxr=1./dxi=1ip1=2im1=1if ((a(ip1,j,k).ne.mdv).and.(a(im1,j,k).ne.mdv)) thend(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))elsed(i,j,k)=mdvendifi=ieip1=ieim1=ie-1if ((a(ip1,j,k).ne.mdv).and.(a(im1,j,k).ne.mdv)) thend(i,j,k)=dxr*(a(ip1,j,k)-a(im1,j,k))elsed(i,j,k)=mdvendifendifc completed gridpoints at i=1, i=ie, j1<=j<=j2enddo ! j-loopc completed gridpoints at 1<j<jec do gridpoints at j=jeif (northpl.eq.1) thenc for these gridpoints, the derivative in the x-direction is ac derivative in the y-direction at another pole-gridpointdy =rerd*dlatdyr=1./(2.*dy)j=jejp1=je-1jm1=je-1do i=1,ieip=mod(i-1+ (ie-1)/4,ie)+1im=mod(i-1+3*(ie-1)/4,ie)+1if ((a(ip,jp1,k).ne.mdv).and.(a(im,jm1,k).ne.mdv))> thend(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))elsed(i,j,k)=mdvendifenddo ! i-loopc completed gridpoints at j=jeendifc do gridpoints at j=1if (southpl.eq.1) thendy =rerd*dlatdyr=1./(2.*dy)j=1jp1=2jm1=2do i=1,ieip=mod(i-1+ (ie-1)/4,ie)+1im=mod(i-1+3*(ie-1)/4,ie)+1if ((a(ip,jp1,k).ne.mdv).and.(a(im,jm1,k).ne.mdv))> thend(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))elsed(i,j,k)=mdvendifenddo ! i-loopendifc completed gridpoints at j=1enddo ! k-loopelse if (dir(1:1).eq.'Y') thenc -----------------------------c derivation in the y-directionc -----------------------------dy =dlat*rerddyr=1./(2.*dy)do k=1,kedo i=1,iec do regular gridpointsdo j=2,je-1jp1=j+1jm1=j-1if ((a(i,jp1,k).ne.mdv).and.(a(i,jm1,k).ne.mdv))> thend(i,j,k)=dyr*(a(i,jp1,k)-a(i,jm1,k))elsed(i,j,k)=mdvendifenddoc do gridpoints at j=jeif (northpl.eq.1) thenc pole-treatmentj=jejm1=j-1jp1=j-1ip=mod(i-1+(ie-1)/2,ie)+1im=iif ((a(ip,jp1,k).ne.mdv).and.(a(im,jm1,k).ne.mdv))> thend(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))elsed(i,j,k)=mdvendifelsec one-sided derivativej=jejm1=j-1jp1=jif ((a(i,jp1,k).ne.mdv).and.(a(i,jm1,k).ne.mdv))> thend(i,j,k)=2.*dyr*(a(i,jp1,k)-a(i,jm1,k))elsed(i,j,k)=mdvendifendifc completed gridpoints at j=jec do gridpoints at j=1if (southpl.eq.1) thenc pole-treatmentj=1jm1=2jp1=2ip=iim=mod(i-1+(ie-1)/2,ie)+1if ((a(ip,jp1,k).ne.mdv).and.(a(im,jm1,k).ne.mdv))> thend(i,j,k)=dyr*(a(ip,jp1,k)-a(im,jm1,k))elsed(i,j,k)=mdvendifelsec one-sided derivativej=1jm1=1jp1=2if ((a(i,jp1,k).ne.mdv).and.(a(i,jm1,k).ne.mdv))> thend(i,j,k)=2.*dyr*(a(i,jp1,k)-a(i,jm1,k))elsed(i,j,k)=mdvendifendifc completed gridpoints at j=1enddoenddoendifendreal function cosd(arg)c-----------------------------------------------------------------------c compute Cos of an argument in Degree instead of Radianc-----------------------------------------------------------------------real,intent(IN) :: argreal cosreal,parameter :: grad2rad=3.1415926/360.cosd=cos(arg*grad2rad)returnendsubroutine pressure(pr,sp,stag3,ie,je,ke,aklev,bklev,aklay,bklay)c =================================================================c argument declarationinteger ie,je,kereal,intent(OUT) :: pr(ie,je,ke)real,intent(IN) :: sp(ie,je),stag3real,intent(IN) :: aklev(ke),bklev(ke),aklay(ke),bklay(ke)c variable declarationinteger i,j,kreal psrfc statement-functions for the computation of pressurereal prlev,prlayinteger isprlev(is)=aklev(is)+bklev(is)*psrfprlay(is)=aklay(is)+bklay(is)*psrfc computation pressuredo i=1,iedo j=1,jepsrf=sp(i,j)do k=1,keif (stag3.eq.0.) thenpr(i,j,k)=prlev(k)elsepr(i,j,k)=prlay(k)endifenddoenddoenddoendsubroutine pres(pr,sp,ie,je,ke,ak,bk)c =====================================c argument declarationinteger ie,je,kereal,intent(OUT) :: pr(ie,je,ke)real,intent(IN) :: sp(ie,je)real,intent(IN) :: ak(ke),bk(ke)c variable declarationinteger i,j,kc computation pressuredo i=1,iedo j=1,jedo k=1,kepr(i,j,k)=ak(k)+bk(k)*sp(i,j)enddoenddoenddoendsubroutine coslat(cl,pollon,pollat,vmin,dx,dy,ie,je)c ====================================================c argument declarationinteger ie,jereal,intent(OUT) :: cl(ie,je)c variable declarationreal pollon,pollat,vmin(3),vmax(3)real lon,lat,rlon,rlatinteger i,jreal phstophreal pidata pi /3.141592654/C Calculate cos(latitude) array and the coriolis parameterif ((pollon.ne.0.).or.(pollat.ne.90.)) thendo j=1,jerlat=vmin(2)+(j-1)*dydo i=1,ierlon=vmin(1)+(i-1)*dxyphys=phstoph(rlat,rlon,pollat,pollon)C if I use sind(lat in deg): troubles at the N-polecl(i,j)=cos(rlat*pi/180.)enddoenddoelsedo j=1,jelat=vmin(2)+(j-1)*dylat=pi*lat/180.do i=1,iecl(i,j)=cos(lat)enddoenddoendifreturnendsubroutine corpar(f,pollon,pollat,vmin,dx,dy,ie,je)c ====================================================c argument declarationinteger ie,jereal,intent(OUT) :: f(ie,je)c variable declarationreal pollon,pollat,vmin(3),vmax(3)real lon,lat,rlon,rlatinteger i,jreal phstophreal pidata pi /3.141592654/C Calculate cos(latitude) array and the coriolis parameterif ((pollon.ne.0.).or.(pollat.ne.90.)) thendo j=1,jerlat=vmin(2)+(j-1)*dydo i=1,ierlon=vmin(1)+(i-1)*dxyphys=phstoph(rlat,rlon,pollat,pollon)C if I use sind(lat in deg): troubles at the N-polelat=2.*pi*yphys/360.f(i,j)=0.000145444*sin(lat)enddoenddoelsedo j=1,jelat=vmin(2)+(j-1)*dylat=pi*lat/180.do i=1,ief(i,j)=0.000145444*sin(lat)enddoenddoendifreturnendsubroutine vint(intfield,field,sp,pmin,pmax,ie,je,ke,ak,bk)C ===========================================================c argument declarationinteger ie,je,kereal intfield(ie,je),field(ie,je,ke),sp(ie,je)real ak(ke),bk(ke)real pmin,pmaxc variable declarationinteger i,j,k,kmin,kmaxreal coeff,pvalc statement-functions for the computation of pressurereal pp,psrfinteger ispp(is)=ak(is)+bk(is)*psrfC ====================================================================C === begin of main part of this subroutine ========================C ====================================================================do i=1,iedo j=1,jeintfield(i,j)=0.psrf=sp(i,j)if (psrf.lt.pmin) thenintfield(i,j)=-999.98999goto 556endifkmin=1kmax=1do k=2,keif ((pp(k).lt.pmax).and.(pp(k-1).gt.pmax)) kmin=kif ((pp(k).lt.pmin).and.(pp(k-1).gt.pmin)) kmax=k-1enddoc check if field is not equal mdvdo k=kmin,kmax+1if (field(i,j,k).eq.-999.98999) thenintfield(i,j)=-999.98999goto 556endifenddoc pmin and pmax are inbetween same pair of layersc interpolate on intermediate value (pmin+pmax)/2.if (kmin.eq.kmax+1) thenpval=(pmin+pmax)/2.coeff=(pval-pp(kmin-1))/(pp(kmin)-pp(kmin-1))intfield(i,j)=(pmax-pmin)*& (field(i,j,kmin)*coeff+field(i,j,kmax)*(1.-coeff))goto 555endifc one layer is inbetween pmin and pmaxif (kmin.eq.kmax) thenif (psrf.lt.pmax) thenintfield(i,j)=field(i,j,kmin)*(psrf-pmin)elseintfield(i,j)=field(i,j,kmin)*(pmax-pmin)endifgoto 555endifc loop for the interior levelsdo k=kmin+1,kmax-1intfield(i,j)=intfield(i,j)+field(i,j,k)*& 0.5*(pp(k-1)-pp(k+1))enddoc special treatment of the bounding levelsif (kmin.eq.1) thenif (psrf.lt.pmax) thenintfield(i,j)=intfield(i,j)+& field(i,j,1)*(0.5*(pp(1)-pp(2))+psrf-pp(1))elseintfield(i,j)=intfield(i,j)+& field(i,j,1)*(0.5*(pp(1)-pp(2))+pmax-pp(1))endifelsecoeff=(pmax-pp(kmin-1))/(pp(kmin)-pp(kmin-1))intfield(i,j)=intfield(i,j)+& field(i,j,kmin)*0.5*(pmax-pp(kmin+1))+& (field(i,j,kmin)*coeff+field(i,j,kmin-1)*(1.-coeff))*& 0.5*(pmax-pp(kmin))endifif (kmax.eq.ke) thenintfield(i,j)=intfield(i,j)+& field(i,j,ke)*0.5*(pp(ke-1)-pp(ke))elsecoeff=(pmin-pp(kmax+1))/(pp(kmax)-pp(kmax+1))intfield(i,j)=intfield(i,j)+& field(i,j,kmax)*0.5*(pp(kmax-1)-pmin)+& (field(i,j,kmax)*coeff+field(i,j,kmax+1)*(1.-coeff))*& 0.5*(pp(kmax)-pmin)endif555 continueC Calculate mean valueif (psrf.lt.pmax) thenintfield(i,j)=intfield(i,j)/(psrf-pmin)elseintfield(i,j)=intfield(i,j)/(pmax-pmin)endif556 continueenddoenddoreturnendsubroutine nvint(intfield,field,sp,pmi,pma,> ie,je,ke,ak,bk,iflag)C ===========================================cc iflag input =0 average (like SR vint); =1 integrationc argument declarationinteger ie,je,ke,iflagreal intfield(ie,je),field(ie,je,ke),sp(ie,je)real ak(ke),bk(ke)real pmin,pmax,pmi,pmac variable declarationinteger ilevinteger i,j,k,kmin,kmax,kshiftreal mdv,coeff,pvalc statement-functions for the computation of pressurereal pp,psrfinteger ispp(is)=ak(is)+bk(is)*psrfC ====================================================================C === begin of main part of this subroutine ========================C ====================================================================C set misdat valuemdv=-999.98999C determine whether data is on model levels (ilev=0)C or pressure levels (ilev=1)ilev=1do k=1,keif (bk(k).ne.0.) ilev=0enddodo i=1,iedo j=1,jeintfield(i,j)=0.psrf=sp(i,j)kmin=1kmax=1pmin=pmipmax=pmaif (ilev.eq.0) thenif (psrf.lt.pmin) thenintfield(i,j)=mdvgoto 456endifendifdo k=2,keif ((pp(k).le.pmax).and.(pp(k-1).gt.pmax)) kmin=kif ((pp(k).le.pmin).and.(pp(k-1).gt.pmin)) kmax=k-1enddoc if model levels: set vint=mdv if one level is equal mdvif (ilev.eq.0) thendo k=kmin,kmax+1if (field(i,j,k).eq.mdv) thenintfield(i,j)=mdvgoto 456endifenddoendifc if pressure levels: check whether lowest levels are not mdvc if they are: adjust pmaxif (ilev.eq.1) thenkshift=0do k=kmin,kmaxif (field(i,j,k).eq.mdv) thenkshift=kshift+1pmax=0.5*(pp(k)+pp(k+1))endifenddoc if mdv occur go even one level higherif (kshift.gt.0) kshift=kshift+1kmin=kmin+kshiftif (kmin.gt.kmax+1) thenintfield(i,j)=mdvgoto 456endifendifc pmin and pmax are inbetween same pair of layersc interpolate on intermediate value (pmin+pmax)/2.if (kmin.eq.kmax+1) thenpval=(pmin+pmax)/2.coeff=(pval-pp(kmin-1))/(pp(kmin)-pp(kmin-1))intfield(i,j)=(pmax-pmin)*& (field(i,j,kmin)*coeff+field(i,j,kmax)*(1.-coeff))goto 455endifc one layer is inbetween pmin and pmaxif (kmin.eq.kmax) thenif (psrf.lt.pmax) thenintfield(i,j)=field(i,j,kmin)*(psrf-pmin)elseintfield(i,j)=field(i,j,kmin)*(pmax-pmin)endifgoto 455endifc loop for the interior levelsdo k=kmin+1,kmax-1intfield(i,j)=intfield(i,j)+field(i,j,k)*& 0.5*(pp(k-1)-pp(k+1))enddoc special treatment of the bounding levelsif (kmin.eq.1) thenif (psrf.lt.pmax) thenintfield(i,j)=intfield(i,j)+& field(i,j,1)*(0.5*(pp(1)-pp(2))+psrf-pp(1))elseintfield(i,j)=intfield(i,j)+& field(i,j,1)*(0.5*(pp(1)-pp(2))+pmax-pp(1))endifelsecoeff=(pmax-pp(kmin-1))/(pp(kmin)-pp(kmin-1))intfield(i,j)=intfield(i,j)+& field(i,j,kmin)*0.5*(pmax-pp(kmin+1))+& (field(i,j,kmin)*coeff+field(i,j,kmin-1)*(1.-coeff))*& 0.5*(pmax-pp(kmin))endifif (kmax.eq.ke) thenintfield(i,j)=intfield(i,j)+& field(i,j,ke)*0.5*(pp(ke-1)-pp(ke))elsecoeff=(pmin-pp(kmax+1))/(pp(kmax)-pp(kmax+1))intfield(i,j)=intfield(i,j)+& field(i,j,kmax)*0.5*(pp(kmax-1)-pmin)+& (field(i,j,kmax)*coeff+field(i,j,kmax+1)*(1.-coeff))*& 0.5*(pp(kmax)-pmin)endif455 continueC Calculate mean or integrated valueif (iflag.eq.0) thenif ((ilev.eq.0).and.(psrf.lt.pmax)) thenintfield(i,j)=intfield(i,j)/(psrf-pmin)elseintfield(i,j)=intfield(i,j)/(pmax-pmin)endifelseintfield(i,j)=100./9.8*intfield(i,j)endif456 continueenddoenddoreturnendC =================================================================subroutine z2hno3(z,nmbl,zval,prof,hno3)C =================================================================C Input is a hno3-profile ("prof") in function of "nmbl" Z-levelsC stored in zval.C Output is the "hno3"-value on this "z"implicit noneinteger nmblC Argument declarationreal zval(nmbl),prof(nmbl),hno3,zC Further variables declarationreal reinteger iC Check if z is inside profile valuesif (z.gt.zval(nmbl)) thenprint*,'*** error: z is over the range of HNO3-profile:'print*,z,zval(nmbl)stopendifif (z.lt.zval(1)) thenprint*,'*** error: z is below the range of HNO3-profile:'print*,z,zval(1)stopendifC Interpolate z of table to find the HNO3 valuedo i=1,nmbl-1if (z.ge.zval(i) .and. z.lt.zval(i+1)) thenre=(z-zval(i))/(zval(i+1)-zval(i))goto 22endifif (i.eq.nmbl-1) thenprint*,'***Error: z:',z,' outside HNO3-Profile...'stopendifenddo22 hno3=prof(i)+re*(prof(i+1)-prof(i))* print*,zi,zval(i),zval(i+1),rereturnendC =================================================================subroutine p2h2o(p,xih2o)C =================================================================implicit noneinteger nnparameter (nn=12)real wprof(nn),pprof(nn)c variables declarationreal p,xih2o,kinteger ndata pprof/150.,130.,100.,90.,80.,70.,60.,50.,40.,30.,20.,10./data wprof/4.57,4.4,4.25,4.28,4.38,4.53,4.7,4.97,5.4,> 5.9,6.6,6.85/do n=1,nn-1if (p.le.pprof(n).and.p.ge.pprof(n+1)) thenk=(p-pprof(n))/(pprof(n+1)-pprof(n))goto 100endifif (n.eq.nn-1) thenprint*,p,pprof(1),pprof(nn)stop'Table too small'endifenddo100 continuexih2o=wprof(n)+k*(wprof(n+1)-wprof(n))returnendC =================================================================subroutine sg(temp,pres,rad,as,xHNO3,xH2O,sedi,growth)C =================================================================C Input of sg are a Temperature[K], a pressure[hPa], a Radius[um],C a particle shape specification stored in "as" with s forC spherical particles, a for aspherical a mixing ratio of HNO3.C Output is the sedimentation factor and a growth rate for theC specified particle.implicit noneC Constants declarationreal deltarho, g, cunn_a, cunn_b, cunn_creal ff1, ff2real gas_diffusivity_corr, mb2mmreal aspect_ratio, d2, k, pi, Rparameter (g=9.81, deltarho=1.62e3) ! deltarho=(NAT Density-Air Density)parameter (cunn_a=1.257, cunn_b=0.4, cunn_c=1.1)!cunn stays for Cunningham...parameter (ff1=1.12, ff2=0.58248) ! Form factors for asphericity 1:3parameter (mb2mm=760./1000.)C For aspher. part., changing the next value implies changing others as formfactorsparameter (aspect_ratio=1./3.)parameter (d2=3.e-10, k=1.380661e-23, pi=3.1415927, R=8.31441)C Arguments declarationdouble precision sedi,growthreal temp,pres,rad,xHNO3,xH2Ocharacter*(1) asC Functions declarationdouble precision cunn_corr,grownC Further variablesdouble precision cc,wStokes,lamda,knudsen,capacity,deltapressreal eta,p_HNO3,p_H2Oeta = 6.45e-8 * tempwStokes = 2.* rad**2 * g * deltarho / (9. * eta)lamda= k * temp / ( 2.**0.5 * 100. * pres * pi * d2**2 )knudsen= lamda / radif (as.eq.'s') thencc=cunn_corr(knudsen,1.,1.,cunn_a,cunn_b,cunn_c)capacity=1.elsecc=cunn_corr(knudsen,ff1,ff2,cunn_a,cunn_b,cunn_c)capacity=sqrt(1.-aspect_ratio**2) / (aspect_ratio**(2./3.) *> log( (1. + sqrt(1.-aspect_ratio**2)) / aspect_ratio) )endifsedi= wStokes * cc / rad**2p_H2O= pres * xH2O * 1.e-6p_HNO3= pres * xHNO3 * 1.e-9deltapress= p_HNO3 - ( exp( (-2.7836 - 0.00088 * temp) *> log (p_H2O * mb2mm) + 89.7674 - 26242. / temp +> 0.021135 * temp ) / mb2mm )growth= capacity * grown(temp,pres,rad,sedi,deltapress,eta)C print*,'T,p,r,shape,xHNO3 -> sedi,growth,',temp,pres,rad,as,C > xHNO3,sedi,growthendC =================================================================double precision function grown(temp,pres,rad,sedi,deltapress,eta)C =================================================================implicit nonereal molmass_g,molmass_s,rho_sreal p0,T0,R,pi,gas_diffusivity_corrparameter (molmass_g=63.e-3 ,molmass_s=117.e-3 ,rho_s=1.62e3 )parameter (p0=1013.25,T0=273.15,pi=3.1415927,R=8.3144)parameter (gas_diffusivity_corr=1.)double precision sedi,deltapress,gasdiffu,gasdiffcorr,prodreal temp,pres,rad,etareal meanvel,Reynold,Schmidt,ventilC if ( temp.lt.233 .or. temp.gt.313 ) print*,'*** Warning:C >Temperature',temp,'is out of validity for diffusivity'gasdiffu= sqrt( 0.018 / molmass_g ) * 0.211e-4 *> ( temp / T0 )**1.94 * (P0 / pres)meanvel= sqrt( 8.* R * temp / (molmass_g * pi) )gasdiffcorr= gasdiffu / ( 1. + 4. * gasdiffu /> ( gas_diffusivity_corr * rad * meanvel ) )Reynold= sedi * rad**2 * 2. * rad / etaSchmidt= eta / gasdiffuprod= Schmidt**(1./3.) * sqrt(Reynold)if ( prod.lt.1.4) thenventil= 1. + 0.108 * prod**2elseventil= 0.78 + 0.308 * prodendifgrown= gasdiffcorr * ventil * molmass_s * deltapress * 100. /> ( rho_s * R * temp )returnendC =================================================================double precision function cunn_corr(knudsen,ff1,ff2,a,b,c)C =================================================================double precision knudsenreal ff1,ff2,a,b,ccunn_corr=(1. + (ff2/ff1) * knudsen * (a+b*exp(-c/knudsen)))/ff1returnend