cxdef 285 linear -25.0 0.33333 cydef 181 linear 0.0 0.33333 ctdef 25 linear 12Z12dec2011 03hr czdef 24 levels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 c 21 22 23 24 ccvars 49 include '/clima/tograds.lib/tograds.inc' character * 8 Cdate8 character * 5 ctime character * 4 clev character * 3 c3,cend,cihr character * 8 ab character * 250 c61 dimension s(8),pradi8(8),pdens8(8),id7(7),ar7(7) real, ALLOCATABLE:: vau (:,:),val (:,:),z(:,:,:),z5(:,:),z7(:,:) real, ALLOCATABLE:: hgt (:,:),dl(:,:,:),sfc(:,:,:),dmsk(:,:) real, ALLOCATABLE:: dd(:,:,:),dw(:,:,:) real, ALLOCATABLE:: s8 (:,:,:,:) real, ALLOCATABLE:: q (:,:,:),qoqi(:,:,:) real, ALLOCATABLE:: pint (:,:,:),sold(:,:,:),zold(:,:,:) real, ALLOCATABLE:: aot (:,:),sfcsvi(:,:),ddsvi(:,:),dwsvi(:,:) real, ALLOCATABLE:: sfc4(:,:),dlsvi(:,:),dl4(:,:) real, ALLOCATABLE:: u10 (:,:),v10(:,:),acp(:,:) real, ALLOCATABLE:: t(:,:,:),tznd(:,:,:),tsrf(:,:,:) + ,fimm(:,:,:),fimm2(:,:,:),dns(:,:,:) + ,fimml(:,:),fimm2l(:,:),psfc(:,:) DATA PRADI8/0.15,0.25,0.47,0.80,1.36,2.29,3.93,7.24/ DATA PDENS8/2.50,2.50,2.50,2.50,2.65,2.65,2.65,2.65/ data id7/13,08,01,0,0,0,0/ wctln(1:1) = ' ' call getarg(1,cdate8) call getarg(2,cend) call getarg(3,cihr) call getarg(4,ab) read (cend,*) iend read (cihr ,*) inchr read (Cdate8(1:2),'(i2)') id7(1) read (Cdate8(3:4),'(i2)') id7(2) id7(3) = 1 call datetohr(id7(3),id7(2),id7(1),id7(4),id7(5),IHRS00) read (Cdate8(5:6),'(i2)') id7(3) read (Cdate8(7:8),'(i2)') id7(4) c print *,' id7=',id7 call datetohr(id7(3),id7(2),id7(1),id7(4),id7(5),IHRS) call RGRADS (id7,81, 001, 0, 1, 1, 1, 0., DDD, NLRET) IMU = rnx_ctl !(1) JMU = rny_ctl !(1) NZ = nz_ctl (1) ALLOCATE(vau(IMU,JMU)) ALLOCATE(val(IMU,JMU)) ALLOCATE(hgt(IMU,JMU)) ALLOCATE(z(IMU,JMU,NZ)) ALLOCATE(q(IMU,JMU,NZ)) ALLOCATE(qoqi(IMU,JMU,NZ)) ALLOCATE(pint(IMU,JMU,NZ)) ALLOCATE(psfc(IMU,JMU)) ALLOCATE(zold(IMU,JMU,NZ+1)) ALLOCATE(sold(IMU,JMU,NZ+1)) ! real, ALLOCATABLE:: t(:,:,:),tznd(:,:,:),tsrf(:,:,:) ! + ,fimm(:,:,:),fimm2(:,:,:),dns(:,:,:) ALLOCATE(t(IMU,JMU,NZ)) ALLOCATE(tznd(IMU,JMU,NZ)) ALLOCATE(tsrf(IMU,JMU,NZ)) ALLOCATE(fimm(IMU,JMU,NZ)) ALLOCATE(fimm2(IMU,JMU,NZ)) ALLOCATE(dns(IMU,JMU,NZ)) ALLOCATE(s8(IMU,JMU,NZ,8)) ALLOCATE(dl(IMU,JMU,8)) ALLOCATE(sfc(IMU,JMU,8)) ALLOCATE(aot(IMU,JMU)) ALLOCATE(dd (IMU,JMU,8)) ALLOCATE(dw (IMU,JMU,8)) ALLOCATE(sfcsvi(IMU,JMU)) ALLOCATE(ddsvi (IMU,JMU)) ALLOCATE(dwsvi (IMU,JMU)) ALLOCATE(dlsvi (IMU,JMU)) ALLOCATE(sfc4(IMU,JMU)) ALLOCATE(dl4(IMU,JMU)) ALLOCATE(u10 (IMU,JMU)) ALLOCATE(v10 (IMU,JMU)) ALLOCATE(acp (IMU,JMU)) ALLOCATE(z5 (IMU,JMU)) ALLOCATE(z7 (IMU,JMU)) ALLOCATE(dmsk(IMU,JMU)) ALLOCATE(fimml(IMU,JMU)) ALLOCATE(fimm2l(IMU,JMU)) ! NOB_CTL(NCTLMX),XOB_CTL(NOBSMX),YOB_CTL(NOBSMX) dlmdU = rxi_ctl !(1) dphdu = ryi_ctl !(1) bwU = rx0_ctl ! (1) bsU = ry0_ctl ! (1) beU = bwU + (IMU-1)*dlmdU bnU = bsU + (JMU-1)*dphdU print *,imu,"f",bwu,beu,dlmdu print *,jmu,"f",bsu,bnu,dphdu print *,NOB_CTL(1) do k=1,NOB_CTL(1) print *,XOB_CTL(k),YOB_CTL(k) enddo c================================= ! open (1,file="AOT.gdat",form="unformatted" c=========================================== open (2,file="xDUST_"//trim(ab)//".gdat" + ,form="unformatted" + ,access="direct",recl=IMU*JMU) open (3,file="xIRON_"//trim(ab)//".gdat" +,form="unformatted" + ,access="direct",recl=IMU*JMU) c================================= do 500 IHT=IHRS+3,IHRS+IEND,inchr c do 500 IHT=IHRS,IHRS+366*24,inchr c do 500 IHT=IHRS,IHRS+31*24,inchr c================================= JHT = IHT call hrtodate(JHT,id7(5),0,id7(3),id7(2),id7(1),id7(4)) call RGRADS (id7,07, 1, 0, IMU, JMU, 1, 0., hgt, NLRET) call RGRADS (id7,07,100,500, IMU, JMU, 1,500., z5, NLRET) call RGRADS (id7,07,100,700, IMU, JMU, 1,700., z7, NLRET) call RGRADS (id7,200,1, 0, IMU, JMU, 1, 0., dmsk, NLRET) call RGRADS (id7,1 ,1, 0, IMU, JMU, 1, 0., psfc, NLRET) do l=1, NZ call RGRADS (id7,7,109,l,IMU,JMU,1,1.*l,z(:,:,l), NLRET) call RGRADS (id7,11,109,l,IMU,JMU,1,1.*l,t(:,:,l), NLRET) call RGRADS (id7,51,109,l,IMU,JMU,1,1.*l,q(:,:,l), NLRET) call RGRADS (id7, 1,109,l,IMU,JMU,1,1.*l,pint(:,:,l), NLRET) do k=1,8 call RGRADS (id7,200+k,109,l,IMU,JMU,1,1.*l,s8(:,:,l,k),NLRET) enddo ! k enddo ! lev cdl idl=221-1 isf=211-1 idd=231-1 idw=241-1 do k=1,8 call RGRADS(id7,idl+k,1,0,IMU,JMU,1,0., dl(:,:,k), NLRET) call RGRADS(id7,isf+k,1,0,IMU,JMU,1,0., sfc(:,:,k), NLRET) call RGRADS(id7,idd+k,1,0,IMU,JMU,1,0., dd (:,:,k), NLRET) call RGRADS(id7,idw+k,1,0,IMU,JMU,1,0., dw (:,:,k), NLRET) enddo call RGRADS(id7,33,105,10,IMU,JMU,1,10.,u10, NLRET) call RGRADS(id7,34,105,10,IMU,JMU,1,10.,v10, NLRET) call RGRADS(id7,61,1,0,IMU,JMU,1,0.,acp, NLRET) cpravi aot! !level 0mnm? do j=1,jmu do i=1,imu do l=1,nz zold(i,j,l) = -hgt(i,j)+z(i,j,l) sold(i,j,l) = 0 do k=1, 8 sold(i,j,l) = sold(i,j,l) +s8(i,j,l,k) enddo enddo enddo enddo do j=1,jmu do i=1,imu dl(i,j,:) = 0 c z(i,j,l) je "ZMID" , u model ZIUNT na interface-u; c dl(i,j,k) = 0 zb = hgt(i,j) do l=1,nz-1 c dz = z(i,j,l)-z0 zt=0.5*(z(i,j,l)+z(i,j,l+1)) dz = zt-zb do k=1,8 dl(i,j,k)=dl(i,j,k)+dz*s8(i,j,l,k) enddo ! k zb = zt enddo enddo ! i enddo ! j aot = 0 sfcsvi = 0 ddsvi = 0 dwsvi = 0 do j=1,jmu do i=1,imu if( hgt(i,j) .lt.5000.and.hgt(i,j).ge.-200 ) then c call checkp(imu,jmu,u10,-100.,100.,und) call dl2aotxx(8, dl(i,j,1:8),aot(i,j)) sfcsvi(i,j) = sum(sfc(i,j,1:8)) dlsvi(i,j) = sum(dl(i,j,1:8)) ddsvi(i,j) = sum(dd(i,j,1:8)) dwsvi(i,j) = sum(dw(i,j,1:8)) sfc4(i,j) = sum(sfc(i,j,1:4)) dl4(i,j) = sum(dl(i,j,1:4)) endif end do end do !iron! do j=1,jmu do i=1,imu ! fimml(i,j) = 0 ! fimm2l(i,j) = 0 ! fimlic(i,j) = 0 ! cwml (i,j) = 0 ! tm55 (i,j) = 0 !--------------------- FIMML(i,j) = 0 FIMM2L(i,j) = 0 zb = hgt(i,j) pb = psfc(i,j) do L=1, NZ-1 zt=0.5*(z(i,j,l)+z(i,j,l+1)) dz = zt-zb !--------------------- TSRF(i,j,l) = 0 TZND(i,j,l) = 0 FIMM(i,j,l) = 0 FIMM2(i,j,l) = 0 DNS (i,j,l) = 0 qoqi (i,j,l) = 0 !--------------------- DO N=1, 7 S(N)=S8(I,J,L,N) ! *1E9 !converted to ug/m^3 ! N particles per unite volume ! ZND(I,J,L,N)=(S(N)*1.E-6) ! */(PDENS8(N)*1.E+6*4./3.*3.14*(PRADI8(N)*1.E-6)**3) !#/m^3 ZND=(S(N)*1.) */(PDENS8(N)*1.E+3*4./3.*3.14*(PRADI8(N)*1.E-6)**3) !#/m^3 TZND(I,J,L)=TZND(I,J,L)+ZND !#/m^3 ! surface concentration for bins SRF=ZND*4*3.14*(PRADI8(N)*1E-6)**2 !m^2/m^3 ! total N of per unite volume ! total surface concentration TSRF(I,J,L)=TSRF(I,J,L)+SRF !m^2/m^3 ENDDO ! n TT=T(I,J,L)-273.2 ! pp=PINT(i,j,l) pp=0.5*(PINT(i,j,l)+pb) CALL tp2qw(pp,tt,qw,qi,qint) qoqi(i,j,l) = q(i,j,l)/qi IF(TT.LT.-5..AND.TT.GT.-55..and.q(i,j,l).ge.0.9*qi)THEN ! IF(TT.LT.-5..AND.TT.GT.-30.)THEN DNS(I,J,L)=EXP(-0.517*TT+8.934) ! [1/m^2] FIMM(I,J,L)=DNS(I,J,L)*TSRF(I,J,L)/1000 ! #/litres ! alternative - modified DeMott (convert tznd to #/cm^3)!FIMM2 in ! #/litres FIMM2(I,J,L)=0.0008*(10**(-0.2*(TT+9.7))) &*(TZND(I,J,L)/1000000.)**1.25 FIMML(i,j) = FIMML(i,j) + 1000.* dz * FIMM(I,J,L) ! #/m2 FIMM2L(i,j) = FIMM2L(i,j) + 1000.*dz * FIMM2(I,J,L) ! #/m2 ! print *,FIMM(I,J,L),i,j,dz endif !--------------------- it=(IHT-IHRS)/inchr+1 !(iyyv-yyst)*12+immv ! if(FIMM2(I,J,L).ne.0) write(77,*),tt,FIMM2(I,J,L) ! ,TT,i,j,l ! if(FIMM2(I,J,L).ne.0) write(78,*),tt,q(i,j,l)/qi ! if(FIMM2(I,J,L).ne.0.and.tt.gt.0) write(7,*),FIMM2(I,J,L),TT,i,j,l zb = zt pb = pint(i,j,l) enddo ! l enddo enddo NZ = 28 ! nlevs NV = 6 !vars NL1 = 8 NL2 = 18 NZ1= (NL2-NL1)+1 jt = (IHT-IHRS00)/inchr do L=NL1, NL2 do iv=1, NV irec=(jt-1)*NZ1*NV+(iv-1)*NZ1+L-8+1 if(iv.eq.1) write(3,rec=irec) fimm(:,:,l) if(iv.eq.2) write(3,rec=irec) fimm2(:,:,l) if(iv.eq.3) write(3,rec=irec) qoqi(:,:,l) if(iv.eq.4) write(3,rec=irec) zold(:,:,l) if(iv.eq.5) write(3,rec=irec) sold(:,:,l) if(iv.eq.6) write(3,rec=irec) t (:,:,l) ! print *,l,iv,irec,"qq" enddo enddo !--------------------- !--------------------- cpisi u aot irec= ((IHT-IHRS)/inchr)*2 cpisi u dust ! irec= ((IHT-IHRS)/inchr)*15 !11 parametra irec= ((IHT-IHRS00)/inchr)*15 !11 parametra und = -999 call checkp(imu,jmu,u10,-100.,100.,und) write(2,rec=irec+1) u10 call checkp(imu,jmu,v10,-100.,100.,und) write(2,rec=irec+2) v10 write(2,rec=irec+3) sfc4 write(2,rec=irec+4) sfcsvi write(2,rec=irec+5) dl4 write(2,rec=irec+6) dlsvi write(2,rec=irec+7) ddsvi write(2,rec=irec+8) dwsvi write(2,rec=irec+9) aot ! call checkp(imu,jmu,acp,-100.,10000.,und) write(2,rec=irec+10) acp ! call checkp(imu,jmu,z5,-100.,10000.,und) write(2,rec=irec+11) z5 ! call checkp(imu,jmu,z7,-100.,10000.,und) write(2,rec=irec+12) z7 write(2,rec=irec+13) dmsk write(2,rec=irec+14) FIMML write(2,rec=irec+15) FIMM2L print *,"fimml=",maxval(FIMML),maxval(FIMML) 500 continue end C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine checkp(im,jm,v,bl,bu,und) dimension v(im,jm) do j=1,jm do i=1,im if( v(i,j).ge.bl.and.v(i,j).lt.bu) then else v(i,j) = und endif enddo enddo return end c================================================================== subroutine dl2aotxx(KPS,dload,aot) dimension QEXT550(8), PRADI(8),PDENS(8),dload(KPS) DATA QEXT550 /1.373,3.303,3.245,2.413,2.262,2.260,2.162,2.108/ DATA PRADI /0.15,0.25,0.45,0.78,1.32,2.24,3.80,7.11/ !microns effective DATA PDENS /2.50,2.50,2.50,2.50,2.65,2.65,2.65,2.65/ !g/cm3 aot8 = 0 do n =1, 8 spart = dload(n) aot8=aot8 + + spart * 1000*3*QEXT550(N)/(4*PRADI(N)*PDENS(N)) enddo aot = aot8 c apradi apdens aqext 2.012500 2.575000 2.390750 c print *,k,1000*spart,OTHIC550 return end !call spline(lold,zslh,vij,y2,lm,zuv,vtil,pp,qq) subroutine spline(nold,xold,yold,y2,nnew,xnew,ynew,p,q) ! ! ****************************************************************** ! * * ! * this is a one-dimensional cubic spline fitting routine * ! * programed for a small scalar machine. * ! * * ! * programer: z. janjic, yugoslav fed. Hydromet. Inst., beograd * ! * * ! * * ! * * ! * nold - number of given values of the function. Must be ge 3. * ! * xold - locations of the points at which the values of the * ! * function are given. Must be in ascending order. * ! * yold - the given values of the function at the points xold. * ! * y2 - the second derivatives at the points xold. If natural * ! * spline is fitted y2(1)=0. And y2(nold)=0. Must be * ! * specified. * ! * nnew - number of values of the function to be calculated. * ! * xnew - locations of the points at which the values of the * ! * function are calculated. Xnew(k) must be ge xold(1) * ! * and le xold(nold). * ! * ynew - the values of the function to be calculated. * ! * p, q - auxiliary vectors of the length nold-2. * ! * * ! ****************************************************************** ! d i m e n s i o n 2 xold(nold),yold(nold),y2(nold),p(nold),q(nold) 3,xnew(nnew),ynew(nnew) !----------------------------------------------------------------------- noldm1=nold-1 ! dxl=xold(2)-xold(1) dxr=xold(3)-xold(2) dydxl=(yold(2)-yold(1))/dxl dydxr=(yold(3)-yold(2))/dxr rtdxc=.5/(dxl+dxr) ! p(1)= rtdxc*(6.*(dydxr-dydxl)-dxl*y2(1)) q(1)=-rtdxc*dxr ! if(nold.eq.3) go to 700 !----------------------------------------------------------------------- k=3 ! 100 dxl=dxr dydxl=dydxr dxr=xold(k+1)-xold(k) dydxr=(yold(k+1)-yold(k))/dxr dxc=dxl+dxr den=1./(dxl*q(k-2)+dxc+dxc) ! p(k-1)= den*(6.*(dydxr-dydxl)-dxl*p(k-2)) q(k-1)=-den*dxr ! k=k+1 if(k.lt.nold) go to 100 !----------------------------------------------------------------------- 700 k=noldm1 ! 200 y2(k)=p(k-1)+q(k-1)*y2(k+1) ! k=k-1 if(k.gt.1) go to 200 !----------------------------------------------------------------------- k1=1 ! 300 xk=xnew(k1) ! do 400 k2=2,nold if(xold(k2).le.xk) go to 400 kold=k2-1 go to 450 400 continue ynew(k1)=yold(nold) go to 600 450 if(k1.eq.1) go to 500 if(k.eq.kold) go to 550 ! 500 k=kold ! y2k=y2(k) y2kp1=y2(k+1) dx=xold(k+1)-xold(k) rdx=1./dx ! ak=.1666667*rdx*(y2kp1-y2k) bk=.5*y2k ck=rdx*(yold(k+1)-yold(k))-.1666667*dx*(y2kp1+y2k+y2k) ! 550 x=xk-xold(k) xsq=x*x ! ynew(k1)=ak*xsq*x+bk*xsq+ck*x+yold(k) ! 600 k1=k1+1 if(k1.le.nnew) go to 300 !----------------------------------------------------------------------- return end subroutine tp2qw(pp,tcel,qw,qi,qint) C*********************************************************************** P A R A M E T E R & (H1=1.E0,H2=2.E0 &, D00=0.E0,D125=.125E0,D5=0.5E0 &, A1=610.78E0,A2=17.2693882E0,A3=273.16E0,A4=35.86E0 &, PQ0=379.90516E0,TRESH=.95E0 &, CP=1004.6E0,ELWV=2.50E6,ELIV=2.834E6,ROW=1.E3,G=9.8E0 &, EPSQ=2.E-12,DLDT=2274.E0,TM10=263.16E0,R=287.0E0 &, CPR=CP*R,RCPR=H1/(CPR)) P A R A M E T E R & (ARCP=A2*(A3-A4)/CP,RCP=H1/CP,PQ0C=PQ0*TRESH,RROG=H1/(ROW*G)) C-------------------- t=273.16+tcel TMT0=T-273.16 TMT15=AMIN1(TMT0,-15.) AI=0.008855 BI=1. IF(TMT0.LT.-20.)THEN AI=0.007225 BI=0.9674 ENDIF QW=PQ0/PP*EXP(A2*(T-A3)/(T-A4)) QI=QW*(BI+AI*AMIN1(TMT0,0.)) QINT=QW*(1.-0.00032*TMT15*(TMT15+15.)) IF(TMT0.LE.-40.)QINT=QI return end