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