      SUBROUTINE TM2(ISTEP,NTSD,IDAT,ICL)
C     ---------------------------------------
C                     14/07/99 
C                     19/04/2000 
C                     04/07/2000 
C                     22/05/2001     
C     ---------------------------------------          
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)      
      PARAMETER (IMJM = IM*JM - JM/2)
      PARAMETER (LP1=LM+1)
      PARAMETER (JAM=6+2*(JM-10))      
      PARAMETER (IMAX=2*IM-1  , JMAX=JM)
      PARAMETER (KKMAX=LM)
      PARAMETER (KMAX=KKMAX)
      PARAMETER (KMIN = LM-KKMAX)
      PARAMETER (KTOP = KMIN + 1)
      PARAMETER (LENIWA=1000)
      REAL PHI0,LAM0

C
      DIMENSION IDAT(3)
C       Date of initial conditions
C     IDAT(1) is day
C     IDAT(2) is month
C     IDAT(3) is year - 1900
C
C-----------------------------------------------------------------------
                             C O M M O N /VRBLS/
     & PD(IMJM),RES(IMJM),FIS(IMJM)
C
     &,U(IMJM,LM),V(IMJM,LM),T(IMJM,LM),Q(IMJM,LM)
C-----------------------------------------------------------------------
                             C O M M O N /MASKS/
     & HBM2  (IMJM),VBM2  (IMJM),VBM3  (IMJM),SM    (IMJM),SICE (IMJM)
C
     &,HTM   (IMJM,LM),VTM   (IMJM,LM)     
C-----------------------------------------------------------------------
                              C O M M O N /MAPOT/
     & LSL,IXM,IYM
     &,ISHDE(49)
     &,R,PT,TSPH
     &,WBD,SBD,TLM0D,TPH0D,DLMD,DPHD,CMLD,DP30,X1P,Y1P
     &,DISLP,Z0SLP,SPL(LSM),ALSL(LSM)
     &,DETA(LM),AETA(LM),ETA(LP1),TSHDE(49)
     &,ERLAM0,CPHI0,SPHI0
     &,KUNIT,LUNIT,MUNIT,IUNIT,NUNIT
C-----------------------------------------------------------------------
                             C O M M O N /LOOPS/
     & KHLA  (JAM),KHHA  (JAM),KVLA  (JAM),KVHA  (JAM)
C
     &,KHL2  (JM),KHH2  (JM),KVL2  (JM),KVH2  (JM)
C
     &,LMH   (IMJM),LMV   (IMJM)     
C------------------------------------------------------------------------ 
      DIMENSION IWA(LENIWA),COBAL(20)
      DIMENSION VERTFL(KKMAX),VERFLX(KKMAX),DFLX(KKMAX+1),
     *          DIFV(IMJM,KMAX)        
      DIMENSION BETA(IMJM,KMAX), IFLAG(IMJM,KMAX) 
      DIMENSION DQDT(IMJM,KMAX)
      DIMENSION COLCON(IMJM), COLAIR(IMJM) 
C      DIMENSION DPREC(IMJM)    
C------------------------------------------------------------------------      
      COMMON/ETATRN/
     * QQL(IMJM,KKMAX),QQH(IMJM,KKMAX),      
     * PS(IMJM,2), QQ(IMJM,KKMAX,2),    
     * UU(IMJM,KKMAX,2), VV(IMJM,KKMAX,2),
     * WW(IMJM,KKMAX,2), RO(IMJM,KKMAX,2),     
     * TT(IMJM,KKMAX,2), ZZ(IMJM,KKMAX+1,2),
     * ALFA(IMJM,KKMAX,2),
     * DIFQ(IMJM,KKMAX),
     * AKH(IMJM,KKMAX,2),AKHS(IMJM,2),
     * TOPG(IMJM)
     *,PREC(IMJM,2),ACPR(IMJM,2)
C
      COMMON/CBT1T2/HINTG(IMJM), HINN(IMJM,KMAX),
     *              HZM(IMJM,KMAX+1),BUFI(IMJM,2),      
     *              QFLUX(IMJM,2) ,EFIREA(IMJM,4),
     *              QSOUR(10),ISOUR(10),JSOUR(10),
     *              XSOUR(10),YSOUR(10),
     *              QFIL(IMJM,KKMAX,2),     
     *              IDATCP(3)
C
C                   IDATCP(3)= IM, ID, IY
C 
      DIMENSION VFLUX(IMJM,2)    
C
      PARAMETER (NPOINT=25)      
C
      COMMON/ETAGRD/
     *      XB,DEX,YB,DEY,
     *      COSETA(JM),XETA(IMAX),YETA(JM),
     *      ITMP(NPOINT), JTMP(NPOINT), KTMP(NPOINT)
      COMMON/ETAREA/CAREA(IMJM)
      COMMON/CITERB/EPS,MAXIT           

      COMMON/TMCTRL/ISTEP1,IBSTEP,DTT,
     *  NSTEP,NSWRI,NSFON,
     *  QMIN,VD,CW,IFDAY,QFLX0,
     *  NSSOUR,QCON0,NSOUR,
     *  IREC,IREC1,IREC2,IREC3,IRECF,
     *  AREAX0,AREAY0,NXAREA,NYAREA,
     *  LOGU,DATE2
      COMMON/CBFIRE/ICLF1,QFLX1,ARFKM2,FLON1,FLAT1,
     *                     FLON2,FLAT2,FLON3,FLAT3         
      COMMON/RECRDS/IREC4                     ! nvp    
          
      CHARACTER*6 DATE2
      CHARACTER*6 DATE3
      LOGICAL LOGU           
      COMMON/TIMECP/t1,t2,t3,t4, tt1, tt2      
      REAL*8 t1,t2,t3,t4, tt1, tt2 
C
      COMMON/TM2LOC/ QKM1, QBAKM1, QAAK, QK, QATK, QAWSK, FQ,
     * DELQ1, DELQ2, DELQ3, DPREC(IMJM)             
C
      DIMENSION WA51(51*51,6), WA51U(51*51,6), WA51V(51*51,6),
     *          WA(IMJM,KKMAX)                      ! nvp 
      COMMON/CFINTY/FINTY(30)        
C-----------------------------------------------------------------------
C      DATA FINTY/0.,0.,0. ,0., 0., 0.,0.,0.,0.,0.,
C     *           0.,0.,0. ,0., 0., 0.,0.,0.,0.,0.,
C     *           0.,0.,0.5,2.0,0.5,0.,0.,0.,0.,0./     
      DATA G /9.81/
      DATA CON/-999.0/ 
C 
      IF(ICL.EQ.0)THEN     
      RETURN
      ELSE      
C     |-------------------------------------------| 
C     |                                           | 
C     |              TIME STEPS START             | 
C     |              ISTEP LOOP START             |
C     |-------------------------------------------|      
      
      RADDEG=3.14159265/180.
      DEG=180./3.14159265                            
      IFILT=0
        
      IBSTEP=0
      call clock(t3)
      
      WRITE(66,*)' ///////////////////////////////////  '
      WRITE(66,*)'#/ ISTEP=',IBSTEP+ISTEP,' # ISTEP=',ISTEP,'/#'
      WRITE(66,*)' ///////////////////////////////////  '
      WRITE(0,*)'#/ ISTEP=',IBSTEP+ISTEP
      
C     ==========================
C          NEXT FIRE AREAS
C     ==========================
      
      IF(MOD(ISTEP,30).EQ.1)THEN
      WRITE(66,*)' FIRE AREAS ON ETA GRID '
      DO K=1,4     
      CALL EXTRAC(EFIREA(1,K),IMJM,AMEAN,AMIN,AMAX,0.,NULL)
      WRITE(66,*)'MEAN,MIN,MAX=',AMEAN,AMIN,AMAX,' NULL=',NULL
      ENDDO
      DO K=1,4
      CALL SUMAR(EFIREA(1,K),IMJM,SUMARI)
      WRITE(66,*)'Total fire area=',SUMARI,' km**2'
      ENDDO        
      ENDIF
Ci       
      IF(MOD(ISTEP,2).NE.0)THEN      
      WRITE(66,*)NAME 
      WRITE(66,*)'NTSD=',NTSD
      WRITE(66,*)IDAT
      ENDIF
                                
      WRITE(66,*)'DIAG. ETADOT ISTEP=',IBSTEP+ISTEP,' */'
      CALL EXTR2M(66,'ETADOT', WW(1,1,2),HTM(1,KTOP),HINTG,
     *          IMJM,KMAX,1)
                  
      CALL AMWEIG(RO(1,1,2),HTM(1,KTOP),PS(1,2),RES,
     *IMJM,KMAX,DETA(KTOP),CAREA,1.) 
             
C |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||                       
C
C                   CALCULATION
C
C |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||  Warning
      DO J=1,IMJM
      DO K=1,KMAX
      BETA(J,K)=0.
      IFLAG(J,K)=0
      ENDDO
      ENDDO
C ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||  Warning

      CALL INMASS(QQ(1,1,1),HTM(1,KTOP),RES,PS(1,1),PT,DETA(KTOP),
     *            KMAX,
     *            TQM1,TAM1,RQM1,TAREA)
      DMQ0=TQM1
      DRQ0=RQM1
C ...........................................................
C
C                    SOURCE CALCULATION
C
C ...........................................................
      
      QFLX=QFLX0                              !  30/11/98
      QFLX=QFLX0*FINTY(MOD(ISTEP-1,30)+1)         !  16/10/01      

      WRITE(66,*)'QFLX=',QFLX,' FOR ISTEP=',ISTEP,
     *' MOD(ISTEP-1,30)=',MOD(ISTEP-1,30),' /#'

      WRITE(66,*)'  =================================='
      WRITE(66,*)'          SURFACE CO2 SOURCES'
      WRITE(66,*)'  =================================='
C                  |=================================|
C                  |                                 |
      WRITE(66,*)'  ======= DISTRIBUTED SOURCES ======'
C                  |                                 |
C                  |=================================|            
      KFIRE=1           
      WRITE(66,*)' DISTRIBUTED CO2 SOURCES ISTEP=',IBSTEP+ISTEP,
     * ' KFIRE=',KFIRE,' #/'
C
C                  ---------------------
C                  FLUX CO2 FROM SURFACE
C                       NEW VERSION
C                  ---------------------
      WRITE(66,*)'---------------------'
      WRITE(66,*)'FLUX CO2 FROM SURFACE'
      WRITE(66,*)'     NEW VERSION'
      WRITE(66,*)'---------------------'
      WRITE(66,*)'IRECF=',IRECF

C =================================================================
C      CALL  RDBUF(4,ISTEP,IRECF,BUFI,EFIREA,IMJM)         
C      CALL SFLUXS(4,ISTEP,VFLUX,QFLX,EFIREA)
C      CALL VFLUXS(4,ISTEP,QFLUX,     EFIREA)      ! ,BUFI,IRECF)
Cx      CALL FLUXON(ISTEP,QFLUX,QFLX,EFIREA)       
C =================================================================
C         From TUPI Version
C         ICLF1 = 0 with one fire only
C         ICLF1 = 1 with distribution fires + one fire
C         ICLF1 = 2 with distribution fires only
C
      QFLX=QFLX0 
           
      IF(ICLF1.GE.1)THEN                                                                        
      CALL RDBUFx(4,ISTEP,IRECF,BUFI,EFIREA,IMJM)                      

C      CALL SFLUXS(4,ISTEP,VFLUX,QFLX,EFIREA)
C      WRITE(66,*)' 1) Loaded mass, ISTEP=',ISTEP,' end/'
C      CALL LOADMS(66,VFLUX,EFIREA,IMJM,DTT)                           

      CALL VFLUXS(4,ISTEP,QFLUX,     EFIREA)      ! ,BUFI,IRECF) 
      WRITE(66,*)'Fires, ISTEP=',ISTEP
      WRITE(66,*)'RDBUF(4,...; SFLUXS(4,...;VFLUXS(4,...'
      CALL EXDIAG(66,'EFIREA(*,1) ',EFIREA,IMJM,1)
      CALL EXDIAG(66,'QFLUX(*,1:2)',QFLUX(1,1),IMJM,2)
      WRITE(66,*)'    '
      WRITE(66,*)' 1) Loaded distribution mass, ISTEP=',ISTEP,' end/'
C23456789 1 23456789 123456789 123456789 123456789 123456789 123456789 1     
      CALL LOADMS(66,QFLUX(1,1),EFIREA,IMJM,DTT)          
      ENDIF     !  ICLF1.GT.1  

      IF(ICLF1.LE.1)THEN 
      WRITE(66,*)'    '                                                  !  16/10/01           
      QFLX2=QFLX1*FINTY(MOD(ISTEP-1,30)+1) 
      CALL NULLI(EFIREA(1,2),IMJM)                ! Clear one fire area
      CALL ONEFIR(EFIREA(1,2),ARFKM2,FLON1,FLAT1) ! add one fire
C only 1      CALL ONEFIR(EFIREA(1,2),ARFKM2,FLON2,FLAT2) ! add one fire                        
      CALL FLUXON(ICLF1,QFLUX,QFLX2,EFIREA(1,2))  ! recalculetes QFLUX
      WRITE(66,*)' 2) Loaded total mass?, ISTEP=',ISTEP,' end/'     
      CALL LOADMS(66,QFLUX(1,1),EFIREA,IMJM,DTT)       
      WRITE(66,*)'FLUXON(ISTEP,...'
      CALL EXDIAG(66,'EFIREA(*,2) ',EFIREA(1,2),IMJM,1)
      CALL EXDIAG(66,'QFLUX(*,1:2)',QFLUX(1,1),IMJM,2)      
      ENDIF     !  ICLF1.LT.3                     
C =================================================================
      WRITE(66,*)'FIRE AREAS ON ETA GRID FOR ISTEP=',ISTEP,' /#'
      DO K=1,4
      CALL SUMAR(EFIREA(1,K),IMJM,SUMARI)
      WRITE(66,*)K,' #Total fire area=',SUMARI,' km**2'
      ENDDO 
      WRITE(66,*)'FIRE AREAS IN CELLS, in km**2'       
      DO K=1,4     
      CALL EXTRAC(EFIREA(1,K),IMJM,AMEAN,AMIN,AMAX,0.,NULL)
      WRITE(66,*)'MEAN=',AMEAN,' MIN=',AMIN,' MAX=',AMAX,' NULL=',NULL
      ENDDO             
C
C ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
C               Distribution sources diagnosis
C ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||      
      CALL EXTRAC(QFLUX,IMJM,AMEAN,AMIN,AMAX,0.,NULL)
      
      CONV1=1.E4*86400.*1.E-6   !  Convert from (m**2 s) to (ha day)

      WRITE(66,*)'AMEAN=',AMEAN*1.E-6,
     *' AMIN=',AMIN*1.E-6,' AMAX=',AMAX*1.E-6,' kg/m**2/s'
      WRITE(66,*)'NULL=',NULL,' #/'
      WRITE(66,*)'AMEAN=',AMEAN*CONV1,
     *' AMIN=',AMIN*CONV1,' AMAX=',AMAX*CONV1,' kg/ha/day'
            
      CALL SCLINP(QFLUX,CAREA,HZM,IMJM,1.E-6,SUMFLX)
      WRITE(66,*)'1) DISTRIBUTED CO2 FLUX FROM SURFACE=',SUMFLX/1000.,
     * ' t/s'
     
      WRITE(66,*)'2) DISTRIBUTED CO2 FLUX FROM SURFACE=',
     * DTT*SUMFLX/1000.0,' t/step  1 step=',DTT,' sec' 
         
      QFLUX0=DTT*QFLX*1.E+6
      CALL SCLINP(EFIREA(1,KFIRE),HZM,HZM,IMJM,QFLUX0,SUMFX1)
      WRITE(66,*)'3) DISTRIBUTED CO2 FLUX FROM SURFACE=',
     * SUMFX1/1000.,' t/step'
     
      CALL SCLINP(EFIREA(1,KFIRE),HZM,HZM,IMJM,1.,SUMFX2)
      WRITE(66,*)'4) FIRE AREA=',SUMFX2,' km**2 ',SUMFX2*1.E6,' m**2'
      
CX      STOP
      
      CALL SFLXD(66,'for VFLUXS',VFLUX,EFIREA,CAREA,HZM,
     *   IMJM,QFLX,DTT)
             
      IF(ISTEP.EQ.1)THEN
      WRITE(66,*)'QFLX              =',QFLX,' kg/m**2/s'
      WRITE(66,*)'1.E+6 x QFLX      =',QFLX*1.E+6,' kg/km**2/s'
      WRITE(66,*)'1.E+6 x DTT x QFLX=',QFLUX0,' kg/km**2/step' 
      ENDIF            
C ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||      
C     I========================================I
C     I                                        I
C     I              POINT SOURCES             I
C     I                                        I
C     I========================================I

CX      GOTO9000
CX
CX9000  CONTINUE

      WRITE(66,*)'POINT SOURCES EXCLUDED'
                      
      CALL SCLINP(QFLUX,CAREA,HZM,IMJM,1.E-6,SUMFLX)
      WRITE(66,*)'TOTAL CO2 FLUX FROM SURFACE=',SUMFLX/1000.,
     * ' t/s'       
      
C --------------------------------------------------------------
C            F I R E  S O U R C E S  E S T I M A T I O N
C --------------------------------------------------------------      
      CALL INMASS(QQ(1,1,1),HTM(1,KTOP),RES,PS(1,1),PT,DETA(KTOP),
     *            KMAX,
     *            TQM2,TAM2,RQM2,TAREA)
      WRITE(66,*)'1) FIRE SOURCES STEP=',IBSTEP+ISTEP,' */'
      WRITE(66,*)' TQM1=',TQM1,' kg x1E6, RQM1=',RQM1,' */'
      WRITE(66,*)' TQM2=',TQM2,' kg x1E6, RQM2=',RQM2,' */'
      WRITE(66,*)'2) FIRE SOURCES STEP=',IBSTEP+ISTEP,' */'
      WRITE(66,*)' TQM2-TQM1=',TQM2-TQM1,' kg x1E6'
      WRITE(66,*)' RQM2-RQM1=',RQM2-RQM1,' */' 
      DMQF=TQM2 
      DRQF=RQM2     
C --------------------------------------------------------------
      QKM1=TQM2              !  kg x 1.E+6  with fire sources 
      
      WRITE(66,*)'AFTER FIRE ESTIMATION QKM1=',QKM1 
      
C      GOTO 7777    
C --------------------------------------------------------------      
                       
      IF(IBSTEP+ISTEP.EQ.1)THEN
C
C     RECORD INITIAL CONDITION
C        
      WRITE(66,*)' *** RECORD INITIAL CONDITION***************' 
      WRITE(66,*)'ISTEP=',ISTEP,' TIME=',0,' hrs'           

      WRITE(66,*)' ////////////// START //////////'
      WRITE(66,*)'           CALL WRHIST(60,... ->  WRH6VR(60,'      
      WRITE(66,*)' ////////////// ///// //////////'  
C ----------------------------------------------------
C
C      
C
C ----------------------------------------------------     

      CALL WRH6VR(60,IREC3,QQ,RO(1,1,2),ZZ,UU,VV,              ! nvp
     * HTM(1,KTOP),HZM,VTM(1,KTOP),LMH,                        ! nvp
     * AETA(KTOP),ETA(KTOP),                                   ! nvp
     * WA51,WA51U,WA51V,WA,KMAX,CON)                           ! nvp
       
      CALL WRSUR(60,IREC3,ACPR,HZM,WA51,CON)                   ! nvp               
      WRITE(66,*)'IREC3=',IREC3                                ! nvp          
             
      ENDIF
      
C =====================================================
C         F L U X   C A L C U L A T I N   FOR t=k-1
C
C         ------------x------- t=k
C                    / 
C         ----------o--------- t=k-1
C
C =====================================================
7777  CONTINUE
      WRITE(66,*)'2) QKM1=',QKM1 
 
      CALL FLUXMS(QQ(1,1,1),HTM(1,KTOP),
     * UU(1,1,2),VV(1,1,1),VTM(1,KTOP),DTT,DELQ1)
     
      WRITE(66,*)'3) QKM1=',QKM1      
     
      CALL FLUXV(QQ(1,1,1),HTM(1,KTOP),WW(1,1,1),DTT*CW,
     * DELQ2,DELQ3,UPAR,TOAR)
      
      QBAKM1=QKM1+DELQ1+DELQ2
      DQOUT=DELQ1+DELQ2
      
      WRITE(66,*)'F L U X   C A L C U L A T I O N   FOR t=k-1'
      WRITE(66,*)'DELQ1                  =',DELQ1,' Later. Flux'
      WRITE(66,*)'DELQ2                  =',DELQ2,' vertic. Flux'
      WRITE(66,*)'DQOUT=DELQ1+DELQ2      =',DQOUT
      WRITE(66,*)'QKM1                   =',QKM1
      WRITE(66,*)'QBAKM1=QKM1+DELQ1+DELQ2=',QBAKM1
                  
C ===================================================== 

      CALL POISTA(66,IBSTEP+ISTEP,'Befor TRPORT QQ(1,1,1)',
     * QQ(1,1,1),HTM(1,KTOP),KMAX)
       
CX      WRITE(66,*)'STOP BEFORE TRPORT(...'         
CX      STOP
               
      call clock(tt1)
C
C     TRANSPORT CALCULATION
C  
      CALL TRPORT(DTT,QMIN,CW)
C      CALL CORREC(QQH,QQL,HTM,IMJM,KMAX,QMIN)
C      
      call clock(tt2)
      WRITE(66,*)'*/ TRANSP CP time=',tt2-tt1,
     *' STEP=',IBSTEP+ISTEP,'/*' 
C
      IF(ISTEP.EQ.2)THEN
      WRITE(66,*)'ISTEP =',ISTEP
      WRITE(66,*)'IBSTEP=',IBSTEP      
      WRITE(66,*)'DTT   =',DTT
      WRITE(66,*)'QMIN  =',QMIN
      WRITE(66,*)'VD    =',VD
      WRITE(66,*)'CW    =',CW
      WRITE(66,*)'LOGU  =',LOGU
      WRITE(66,*)'QBAKM1=',QBAKM1
      
      CALL EXTRAC(HINTG,IMJM,AMEAN,AMIN,AMAX,0.,NULL)
      WRITE(66,*)'HINTG:  AMEAN,AMIN,AMAX,NULL=',
     * AMEAN,AMIN,AMAX,NULL 
      CALL EXTRAC(HINN,IMJM*KMAX,AMEAN,AMIN,AMAX,0.,NULL)
      WRITE(66,*)'HINN:  AMEAN,AMIN,AMAX,NULL=',
     * AMEAN,AMIN,AMAX,NULL
      CALL EXTRAC(HTM,IMJM*KMAX,AMEAN,AMIN,AMAX,0.,NULL)
      WRITE(66,*)'HTM:  AMEAN,AMIN,AMAX,NULL=',
     * AMEAN,AMIN,AMAX,NULL              

C      STOP
      
      ENDIF
C
      CALL POISTA(66,IBSTEP+ISTEP,'After TRANSP, QQ(1,1,2)',
     * QQ(1,1,2),HTM(1,KTOP),KMAX) 
     
      CALL INMASS(QQ(1,1,2),HTM(1,KTOP),RES,PS(1,2),PT,
     *     DETA(KTOP),KMAX,
     *     DMQA,TAM2,DRQA,TAREA)
C ================================== 
      QAAK=DMQA       ! After advection
C =====================================================
C         F L U X   C A L C U L A T I N   FOR t=k
C
C         ------------x------- t=k
C                    / 
C         ----------o--------- t=k-1
C
C =====================================================
 
      CALL FLUXMS(QQ(1,1,2),HTM(1,KTOP),
     * UU(1,1,2),VV(1,1,2),VTM(1,KTOP),DTT,DELQ1)
     
      CALL FLUXV(QQ(1,1,2),HTM(1,KTOP),WW(1,1,2),DTT*CW,
     * DELQ2,DELQ3,UPAR,TOAR)
     
CX      QK=QKM1+DELQ1+DELQ2
      
      QK=QKM1+0.5*(DELQ1+DELQ2+DQOUT)       !  22/05/2001
      QBAKM1=QKM1+0.5*(DELQ1+DELQ2+DQOUT)   !  22/05/2001         
C =====================================================
C                                                  From To
C                                                    |   |
      CALL VLFLUX(QQ(1,1,2),HTM(1,KTOP),WW(1,1,2),CW,1,KMAX,DTT,
     * VERTFL,VERFLX)
      WRITE(66,*)'VERTICAL FLUX1 STEP= ',IBSTEP+ISTEP,' */'
      DO KL=1,KMAX-1
      WRITE(66,*)KL,VERTFL(KL),' 1.E+6 kg/step: ^ plus; V minus'
      ENDDO
C =====================================================      
      WRITE(66,*)'FLUX STEP=',IBSTEP+ISTEP,
     * '/* HOR=',DELQ1,' VER=',DELQ2,'*/' 
      WRITE(66,*)'FLUX = a)',DELQ2,'b)',DELQ3,
     *'c)',UPAR/1.E6,'d)',TOAR/1.E6,'*/'
      WRITE(66,*)'OUTFLUX STEP=',IBSTEP+ISTEP,
     * ' DQOUT=',DQOUT,' */'
      WRITE(66,*)'DELQ1 - the lateral flux  =',DELQ1
      WRITE(66,*)'DELQ2 - the vertical flux =',DELQ2
C
C      DELQ1 - the lateral flux    
C      DELQ2 - the vertical flux   kg x 1.E+6
C      DELQ3 - the normalazed vertical flux = total flux / area 
C                                 (kg x 1.E+6)/m**2
C       UPAR - the inner area with up velocity in km**2
C       TOAR - the total inner area in km**2
CCCCCC ================================================= CCCCC
CCCCCC         CONSERVATION WITHOUT V.TURBULENCE         CCCCC
CCCCCC ================================================= CCCCC
C     ----------------------------
C     QUASI-CONSERVATION ALGORITHM
C     ----------------------------
      WRITE(66,*)'QUASI-CONSERVATION ALGORITHM'

C      STOP

      call clock(tt1)     
      FQ=QBAKM1
      
      WRITE(66,*)'FQ=QBAKM1 =',FQ,' */'
      
      CALL CONSV1(QQH,QQL,RO,HINN,IMJM*KMAX,
     *     ALFA,BETA,FQ,CSTAR,IR)

      IF(IR.EQ.0)THEN
      ICON=1
      CALL CONSV2(QQH,QQL,QQ(1,1,2),IMJM*KMAX,
     *      ALFA,BETA,IFLAG,CSTAR,ICON)

      IF(ICON.EQ.0)THEN
      CALL CONMAS(QQ(1,1,2),HTM(1,KTOP),RES,PS(1,2),PT,
     * DETA(KTOP),KMAX,0.,QMIN,
     * FQ,DQADJ,QMNEW)
      ENDIF      
   
      ELSE
      CALL CONSV3(QQH,QQL,QQ(1,1,2),IMJM*KMAX,
     *     ALFA)
      ENDIF      
      
      call clock(tt2)
      WRITE(66,*)'*/ CONSERV CP time=',tt2-tt1,
     *' STEP=',IBSTEP+ISTEP,'/*'

CX      WRITE(66,*)'STOP AFTER CONSV'
CX      STOP
      
      CALL IPOSTA(66,IBSTEP+ISTEP,'BETA',BETA,HTM(1,KTOP),KMAX)
      CALL IPOSTA(66,IBSTEP+ISTEP,'ALFA',ALFA,HTM(1,KTOP),KMAX)
         WRITE(66,*)' ****** CONSERVATION, ISTEP=',IBSTEP+ISTEP
         WRITE(66,*)'FQ   =',FQ
         WRITE(66,*)'CSTAR=',CSTAR
         WRITE(66,*)'IR   =',IR      
      
      CALL POISTA(66,IBSTEP+ISTEP,'After CONSERV., QQ(1,1,2)',
     * QQ(1,1,2),HTM(1,KTOP),KMAX)

      CALL INMASS(QQ(1,1,2),HTM(1,KTOP),RES,PS(1,2),PT,
     *            DETA(KTOP),KMAX,
     *            TQM,TAM,RQM,TAREA)

      QACK=TQM      
      DMQT=TQM
      DRQT=RQM
       
      WRITE(66,*)'CONSERVATION=============='
      WRITE(66,*)'AFTER ADVECTION=',QBAKM1,' */'
      WRITE(66,*)'AFTER CONSERV  =',TQM,' */'   
      
CX      WRITE(66,*)'STOP BEFORE VERTICAL DIFFUSION'
CX      STOP
        
CCCCCC ================================================= CCCCC 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC                             
C
C     VERTICAL DIFFUSION
C      
      call clock(tt1)
C ---------------------------------------------------
   
      CALL QDIFF(QQ(1,1,2),QFLUX,DFLX,DIFV,LMH,AKH,
     * RO(1,1,2),ZZ(1,1,2),KMAX,DTT)
      
      CALL INMASS(QQ(1,1,2),HTM(1,KTOP),RES,PS(1,2),PT,
     *            DETA(KTOP),KMAX,
     *            TQM,TAM,RQM,TAREA)       
      QATK=TQM      
      DMQT=TQM
      DRQT=RQM    
C ---------------------------------------------------            
      call clock(tt2)
      WRITE(66,*)'*/ VSDIFF CP time=',tt2-tt1,
     *' STEP=',IBSTEP+ISTEP,'/*' 
           
      WRITE(66,*)'VERTICAL FLUX2 STEP= ',IBSTEP+ISTEP,' */'
      DO KL=1,KMAX+1
      WRITE(66,*)KL,DFLX(KL),' 1.E+6 kg/step: ^ plus; V minus'
      ENDDO 
          
      CALL POISTA(66,IBSTEP+ISTEP,'After VSDIFF, QQ(1,1,2)',
     * QQ(1,1,2),HTM(1,KTOP),KMAX)

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C                      WET SEDIMENTATION
C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
      WRITE(0,*)'NO WET SEDIMENTATION'
      WRITE(0,*)'ISTEP=',ISTEP 
     
      CALL SUBTRA(ACPR(1,2),ACPR(1,1),IMJM,DPREC)
      CALL EXTRAC(DPREC,IMJM,AMEAN,AMIN,AMAX,0.,NULL)
      WRITE(66,*)'DPREC: AMEAN,AMIN,AMAX,NULL=',
     * AMEAN,AMIN,AMAX,NULL       
      
      WRITE(66,*)'SUBTRA... O.K.'
      COEWS=2.8
      COEWS=1.4            
CNO      CALL WTSCA(QQ(1,1,2),COEWS,ACPR,DQDT,HTM(1,KTOP),KMAX,DTT)
C      CALL WTSCA(QQ(1,1,2),2.8,ACPR,DQDT,HTM(1,KTOP),KMAX,DTT)
C                           0.1      
                         
C      CALL WETSED(QQ(1,1,2),ACPR,DQDT,HTM(1,KTOP),KMAX,DTT)      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CALL INMASS(QQ(1,1,2),HTM(1,KTOP),RES,PS(1,2),PT,
     *            DETA(KTOP),KMAX,
     *            TQM,TAM,RQM,TAREA)

      QAWSK=TQM      
      DMQT=TQM
      DRQT=RQM
      
      WRITE(66,*)'   MASS BALANCE FOR STEP=',IBSTEP+ISTEP,' */'
      WRITE(66,*)'DMQ0  =',DMQ0  ,' 1.E6 kg '
      WRITE(66,*)'QKM1  =',QKM1  ,' 1.E6 kg '
      WRITE(66,*)'QBAKM1=',QBAKM1,' 1.E6 kg '
      WRITE(66,*)'QAAK  =',QAAK  ,' 1.E6 kg '
      WRITE(66,*)'QACK  =',QACK  ,' 1.E6 kg '
      WRITE(66,*)'QATK  =',QATK  ,' 1.E6 kg ' 
      WRITE(66,*)'QAWSK =',QAWSK ,' 1.E6 kg '      
      WRITE(0,*)'START CLMDOM'       
C
C     * DETA(KTOP),KMAX,
C
CX      WRITE(0,*)'FINISH COLMAS'

      WRITE(66,*)'START COLMAS ISTEP=',ISTEP,' */' 
          
            CALL CLMDOM(QQ(1,1,2),HTM(1,KTOP),RES,PS(1,2),
     *      DETA(KTOP),KMAX,
     *      COLCON,COLAIR)
     
      CALL EXTRAC(COLCON,IMJM,AMEAN,AMIN,AMAX,0.,NULL)
      WRITE(66,*)'COLCON: AMEAN, AMIN, AMAX, NULL=',
     * 1000.*AMEAN,1000.*AMIN,1000.*AMAX,NULL,' g/m^2' 
           
      WRITE(66,*)'FINISH COLMAS'   
     
C19/02/99      CALL COLDIA(6,IBSTEP+ISTEP,
C19/02/99     * COLCON,COLAIR,1.E-6,1.)    ! in kg 


C      STOP
      WRITE(66,*)'AREAS BALANCE CALCULATION JUMP'
      GOTO 777
C     ..........................................
C
C             AREAS BALANCE CALCULATION
C
C     ..........................................
      AREAX0=-66.0
      AREAY0=-14.0      
      CALL PARARE(AREAX0,AREAY0,I0,J0,TLM0D*RADDEG,TPH0D*RADDEG)
            
      WRITE(66,*)'PARARE: I0=',I0,' J0=',J0

      NXAREA=25
      NYAREA=25
      
C/      NXAREA=35
C/      NYAREA=45      
      
C/      NXAREA=5
C/      NYAREA=5

      IF(ISTEP.EQ.1)THEN
      CALL PASSAB(QQ(1,1,1),IMJM*KMAX ,QFIL(1,1,1))
      CALL ENFILT(QQ(1,1,1),QFIL(1,1,1)) 
      ENDIF
      
      
      CALL PASSAB(QQ(1,1,2),IMJM*KMAX ,QFIL(1,1,2)) 
      CALL ENFILT(QQ(1,1,2),QFIL(1,1,2)) 

      WRITE(66,*)'Longitude:',DEX*(NXAREA-1),' degr.'
      WRITE(66,*)'From ',AREAX0,' to ',AREAX0+DEX*(NXAREA-1)
      WRITE(66,*)'Latitude : ',DEY*(NYAREA-1),' degr.'
      WRITE(66,*)'From ',AREAY0,' to ',AREAY0+DEY*(NYAREA-1)      
      
      CALL BALARE(I0,NXAREA,J0,NYAREA,DTT,
     * QFIL(1,1,1),QFLUX,HTM,UU(1,1,1),VV(1,1,1),VTM,PS(1,1),
     * IWA,LENIWA,COBAL)
     
C/      CALL LOCADR(6,IWA,NXAREA,NYAREA)
      WRITE(66,*)'     BALANCE FOR FILTERING DATA'      
      WRITE(66,*)'AREA BALANCE: I=',I0,' J=',J0,' ISTEP=',ISTEP,' #/'
      WRITE(66,*)'              M=',NXAREA,' N=',NYAREA,' #/'
      WRITE(66,*)'1) TOTM=',COBAL(1),' #/ 1.E6 kg'
      WRITE(66,*)'1) FBOT=',COBAL(2),' #/ 1.E6 kg/step'
      WRITE(66,*)'1) FLAT=',COBAL(3),' #/ 1.E6 kg/step'
      WRITE(66,*)'      '
      WRITE(66,*)'1) FSOU=',COBAL(4),' #/ 1.E6 kg/step'
      WRITE(66,*)'1) FEST=',COBAL(5),' #/ 1.E6 kg/step'
      WRITE(66,*)'1) FNOR=',COBAL(6),' #/ 1.E6 kg/step'
      WRITE(66,*)'1) FWES=',COBAL(7),' #/ 1.E6 kg/step' 
      WRITE(66,*)'1) AREA=',COBAL(8),' #/ 1.E6 kg/step'
      WRITE(66,*)'1) QMEAN=',COBAL(9),' #/ ppm'
      WRITE(66,*)'1) COVQ=',COBAL(10),' #/ ppm'
      WRITE(66,*)'1) AIMA=',COBAL(11),' #/ kg'           
      WRITE(66,*)'      '
      WRITE(66,*)'      '            
      CALL BALARE(I0,NXAREA,J0,NYAREA,DTT,
     * QFIL(1,1,2),QFLUX,HTM,UU(1,1,2),VV(1,1,2),VTM,PS(1,2),
     * IWA,LENIWA,COBAL)
      WRITE(66,*)'2) TOTM=',COBAL(1),' #/ 1.E6 kg'
      WRITE(66,*)'2) FBOT=',COBAL(2),' #/ 1.E6 kg/step'
      WRITE(66,*)'2) FLAT=',COBAL(3),' #/ 1.E6 kg/step'
      WRITE(66,*)'      '
      WRITE(66,*)'2) FSOU=',COBAL(4),' #/ 1.E6 kg/step'
      WRITE(66,*)'2) FEST=',COBAL(5),' #/ 1.E6 kg/step'
      WRITE(66,*)'2) FNOR=',COBAL(6),' #/ 1.E6 kg/step'
      WRITE(66,*)'2) FWES=',COBAL(7),' #/ 1.E6 kg/step'                       
      WRITE(66,*)'2) AREA=',COBAL(8),' #/ 1.E6 kg/step'
      WRITE(66,*)'2) QMEAN=',COBAL(9),' #/ ppm'
      WRITE(66,*)'2) COVQ=',COBAL(10),' #/ ppm'         
      WRITE(66,*)'2) AIMA=',COBAL(11),' #/ kg'
                  
      CALL LATFLX(QFIL,HTM,UU,VV,VTM,PS,IWA,LEN,FLXL)
      WRITE(66,*)'3) FLAT ISTEP=',ISTEP,' FLXL=',FLXL*DTT,
     * ' #/ 1.E6 kg/step'      
          
      IF(ISTEP.EQ.1)THEN
      OPEN(50,FILE=
     * '/pesqout/sdg04/igor/Data/dataout/bnd1area.dat',
     *  STATUS='UNKNOWN')       
      CALL BNDAR(50,I0,NXAREA,J0,NYAREA,3)
      WRITE(66,*)'OPEN AND READ of bnd1area.dat done' 
      ENDIF
      
      STOP 
           
777   CONTINUE

      CALL ATMOSM(66,PS,ATM1,ATM2)
     
      CALL PASSAB(QQ(1,1,2),IMJM*KMAX ,QQ(1,1,1))      
      CALL PASSAB(PS(1,2),IMJM ,PS(1,1))   
      CALL PASSAB(TT(1,1,2),IMJM*KMAX ,TT(1,1,1))
      CALL PASSAB(UU(1,1,2),IMJM*KMAX ,UU(1,1,1))                   
      CALL PASSAB(VV(1,1,2),IMJM*KMAX ,VV(1,1,1))
      CALL PASSAB(WW(1,1,2),IMJM*KMAX ,WW(1,1,1))
      CALL PASSAB(RO(1,1,2),IMJM*KMAX ,RO(1,1,1))
      
      CALL PASSAB(ZZ(1,1,2),IMJM*(KMAX+1) ,ZZ(1,1,1))
      CALL PASSAB(AKH(1,1,2),IMJM*KMAX ,AKH(1,1,1))
      CALL PASSAB(AKHS(1,2),IMJM ,AKHS(1,1))
      CALL PASSAB(PREC(1,2),IMJM ,PREC(1,1))
      CALL PASSAB(ACPR(1,2),IMJM ,ACPR(1,1)) 
                                
      CALL PASSAB(QFIL(1,1,2),IMJM*KMAX ,QFIL(1,1,2))
             
      CALL EXTR2M(66,'CO2', QQ(1,1,1),HTM(1,KTOP),HINTG,
     *          IMJM,KMAX,1)
      CALL EXTR2M(66,'VDIFF', DIFV,HTM(1,KTOP),HINTG,
     *          IMJM,KMAX,1)

      CALL EXTR2M(66,'AKH', AKH,HTM(1,KTOP),HINTG,
     *          IMJM,KMAX,1)                 
C     --------------------------
C
C     --------------------------
      CALL AVSUR(ACPR,HTM1,CON,ACPI,ACPA)
      WRITE(66,*)'#/ ACPR ISTEP=',IBSTEP+ISTEP,
     *  'ACPI=',ACPI,' ACPA=',ACPA,' /#  Old'
           
      CALL AVSUR(ACPR,HINTG,CON,ACPI,ACPA)            
      WRITE(66,*)'#/ ACPR ISTEP=',IBSTEP+ISTEP,
     *  'ACPI=',ACPI,' ACPA=',ACPA,' /#'
     
          
      CALL AVSUR(PREC,HTM1,CON,PREI,PREA)
      WRITE(66,*)'#/ PREC ISTEP=',IBSTEP+ISTEP,
     *  'PREI=',PREI,' PREA=',PREA,' /#  Old'

      CALL AVSUR(PREC,HINTG,CON,PREI,PREA)            
      WRITE(66,*)'#/ PREC ISTEP=',IBSTEP+ISTEP,
     *  'PREI=',PREI,' PREA=',PREA,' /#'       
                
      IF(ISTEP.EQ.2)CALL SITIOS(66,WA)
C                                                  !     22/05/2001
      IF(MOD(ISTEP,NSWRI).EQ.0.OR.ISTEP.EQ.53)THEN ! for 42h or 24+18h
C     -----------------------------
C
C       OPTICAL PARAMETER RECORD
C 
C     -----------------------------
      WRITE(66,*)' OPTICAL PARAMETER RECORD, ISTEP=',
     * ISTEP,' #/'  
      CALL SITIOS(66,WA)
C23456
C     -----------------------------
C
C             OUTPUT RECORD
C 
C     ----------------------------- 
      WRITE(66,*)'OUTPUT RECORD ISTEP=',ISTEP  
           
      CALL EXTRAC(DPREC,IMJM,AMEAN,AMIN,AMAX,0.,NULL)
      WRITE(66,*)'DPREC: AMEAN,AMIN,AMAX,NULL=',
     * AMEAN,AMIN,AMAX,NULL
     
      CALL WRH6VR(60,IREC3,QQ,RO,ZZ,UU,VV,            !  nvp
     * HTM(1,KTOP),HZM,VTM(1,KTOP),LMH,               !  nvp
     * AETA(KTOP),ETA(KTOP),                          !  nvp
     * WA51,WA51U,WA51V,WA,KMAX,CON)                  !  nvp
     
            
      WRITE(66,*)'CALL CALL WRH6VR(60,IREC3,QQ,RO,ZZ,UU,VV, ,...+...'
      WRITE(66,*)'IREC3=',IREC3                       !  nvp
                
      CALL WRSUR(60,IREC3,DPREC,HZM,WA51,CON)         !  nvp       
      WRITE(66,*)'IREC3=',IREC3                       !  nvp     
            
      WRITE(66,*)'END OF OUTPUT RECORD ISTEP=',ISTEP      
      ENDIF               
                                            
      call clock(t4)
      write(66,103) t4-t3 
           
      IF(ISTEP.EQ.NSTEP)THEN 
C     -----------------------------
C             CHECK POINT
C     -----------------------------      
      CALL CHKPT(58,'W',IDAT,IBSTEP+ISTEP,
     * IREC,IREC1,IREC2,QQ,IMJM,KMAX)
C23456
      WRITE(66,*)'CHKPT(58,''W'',IDAT,..., write'
    
      call clock(t2)
      write(66,102) t2-t1     
102   format('CPU TIME for calculation: ',f10.4,' sec')
103   format('CPU TIME for step calculation: ',f10.4,' sec') 
             
C      STOP 'ETA-TRANSPORT MODEL'      

      ENDIF      

      RETURN
      
      ENDIF
C     See top
      END 
C& CORREC
      SUBROUTINE CORREC(A,B,H,IMJM,KMAX,QMIN)
      DIMENSION A(IMJM,KMAX),B(IMJM,KMAX),H(IMJM,KMAX)
      
      DO K=1,KMAX
      DO J=1,IMJM
      IF(H(J).EQ.1.)THEN
        IF(A(J)*H(J).LT.QMIN)A(J)=0.
        IF(B(J)*H(J).LT.QMIN)B(J)=0.
      ENDIF
      ENDDO
      ENDDO
      RETURN
      END
C& QDIFF
      SUBROUTINE QDIFF(Q,FSUR,DFLX,DIFLV,LMH,AKH,
     * RO,ZINT,KMAX,DTQ2)
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)      
      PARAMETER (IMJM = IM*JM - JM/2)
      PARAMETER (LP1=LM+1, LM1=LM-1)
      DIMENSION Q(IMJM,KMAX),FSUR(IMJM),LMH(IMJM),
     * AKH(IMJM,KMAX),RO(IMJM,KMAX),ZINT(IMJM,KMAX+1), 
     * DFLX(KMAX+1),DIFLV(IMJM,KMAX)      
C-----------------------------------------------------------------------
      DIMENSION
     & DC(LM1),AC(LM),BC(LM),CC(LM),FC(LM),
     & ROC(LM),ZC(LP1),CKH(LM1),
     & SOLC(LM),BET(LM+1),Z(LM+1),
     & TFLX(LM+1)
      DIMENSION DIRO(LM),DIZ(LP1),DIKH(LM1)
C                             C O M M O N
      COMMON/INNERP/IPN,INP(IMJM),IINP(IMJM),JINP(IMJM)      
      COMMON/ETAREA/CAREA(IMJM)
C
C                     S T A R T
C

CX      WRITE(66,*)'THE QDIFF SUBROUTINE STARTS'
CX      STOP
      
      DO K=1,LM+1
      TFLX(K)=0.
      ENDDO
                  
      DO1000 J=1,IPN
      JP=INP(J)
      KB=LMH(JP)
      DS=CAREA(JP)
            
CX         IF(J.GT.22680)THEN
CX         WRITE(0,*)'J=',J,' JP=',JP,' KB=',KB 
CX         ENDIF
                     
      ZC(1)=ZINT(JP,1)
      DO K=1,KB
      ZC(K+1)=ZINT(JP,K+1)
      ROC(K)=RO(JP,K)/DS/(ZC(K)-ZC(K+1))   !   21/02/99      
      ENDDO
      
      DO K=1,KB-1
      CKH(K)=AKH(JP,K)
      ENDDO
      
      DO K=1,KB-1
      DC(K)=(ROC(K)+ROC(K+1))*CKH(K)/(ZC(K)-ZC(K+2))
      FC(K)=Q(JP,K)
      ENDDO
      FC(KB)=Q(JP,KB)+DTQ2*FSUR(JP)/(ROC(KB)*(ZC(KB)-ZC(KB+1)))
C                    -
      
      AC(1)=0.
      DO K=2,KB          !   -1
      AC(K)=DTQ2*DC(K-1)/(ROC(K)*(ZC(K)-ZC(K+1)))
      ENDDO
      
      DO K=1,KB-1
      CC(K)=DTQ2*DC(K)/(ROC(K)*(ZC(K)-ZC(K+1)))
      ENDDO
      CC(KB)=0.
      
      DO K=1,KB
      BC(K)=1.+AC(K)+CC(K)
      ENDDO

CX      WRITE(66,*)'STOP BEFORE PROGON'
CX      STOP
CX      IF(J.EQ.5590)THEN
CX      WRITE(66,*)' //////////////////////'
CX      WRITE(66,*)'   POINT=',J
CX      WRITE(66,*)' //////////////////////'
CX      WRITE(66,*)'FSUR(JP)=',FSUR(JP)
CX      WRITE(66,*)' K    q    Ro    z    KH'   
CX      DO KL=1,KB
CX      WRITE(66,*)KL,Q(JP,KL),ROC(KL),ZC(KL),CKH(KL)
CX      ENDDO
CX      WRITE(66,*)'ZC(KB+1)=',ZC(KB+1)
CX                  
CX      ENDIF
      
      CALL PROGON(AC,BC,CC,FC,KB,SOLC,BET,Z)
C23456      
CX      IF(J.EQ.5590)THEN 
CX      WRITE(66,*)'      SOLUTION'     
CX      DO KL=1,KB      
CX      WRITE(66,*)KL,SOLC(KL)
CX      ENDDO
CX      
CX      WRITE(66,*)'         EQUATION'
CX      WRITE(66,*)'KL,AC(KL),BC(KL),CC(KL),FC(KL),DC(KL)'     
CX      DO KL=1,KB      
CX      WRITE(66,*)KL,AC(KL),BC(KL),CC(KL),FC(KL),DC(KL)
CX      ENDDO       
CX
CX      ENDIF
            
      DO K=1,KB
      DIFLV(JP,K)=SOLC(K)-Q(JP,K)
      Q(JP,K)=SOLC(K)
      ENDDO
C
C             TOTAL FLUX OVER ETA SURFACES
C 
      TFLX(1)=0.    
      DO K=1,KB-1
      TFLX(K+1)=TFLX(K+1)+DS*DC(K)*(SOLC(K)-SOLC(K+1))
      ENDDO
      TFLX(KB+1)=TFLX(KB+1)+DS*FSUR(JP)
C                          -      
      
      
1000  CONTINUE                            !     J loop

CX      WRITE(0,*)'End of 1000 loop'
C     Let us calculate the turbulent flux !

      DO K=1,LM+1
      DFLX(K)=DTQ2*TFLX(K)
      ENDDO
      
      RETURN
      END            
C *********************************************
C                                             *
C            FROM    TRPORT.f                 *
C                                             *
C ********************************************* 
C
C     Version 1/12/98
C
C& TRPORT
      SUBROUTINE TRPORT(DT,QMIN,CW)
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)
      PARAMETER (LDA=LM+9,LA=13)               
C
      PARAMETER (
     &  IMJM=IM*JM-JM/2,JMP1=JM+1,JAM=6+2*(JM-10),LB=2*IM+JM-3
     &, LM1=LM-1,LP1=LM+1,L3=LDA+LM+1,L4=LDA+2*LM+1,L7=L4+3
     &, LSCRCH=4*LM+1+LA+1)
     
      PARAMETER (IMT=2*IM-1)      
      PARAMETER (IMAX=2*IM-1  , JMAX=JM)
      PARAMETER (KKMAX=LM)
      PARAMETER (KMIN =LM-KKMAX)
      PARAMETER (KMAX=KKMAX)
      PARAMETER (KTOP=KMIN+1)

      COMMON/ETATRN/
     * QQL(IMJM,KKMAX),QQH(IMJM,KKMAX),      
     * PS(IMJM,2), QQ(IMJM,KKMAX,2),   
     * UU(IMJM,KKMAX,2), VV(IMJM,KKMAX,2),
     * WW(IMJM,KKMAX,2), RO(IMJM,KKMAX,2),     
     * TT(IMJM,KKMAX,2), ZZ(IMJM,KKMAX+1,2),
     * ALFA(IMJM,KKMAX,2),
     * DIFQ(IMJM,KKMAX),
     * AKH(IMJM,KKMAX,2),AKHS(IMJM,2),
     * TOPG(IMJM)
     *,PREC(IMJM,2),ACPR(IMJM,2)
     
      COMMON/DBGXYZ/IXYZ 
                      
      PARAMETER (NPOINT=25)      
C
      COMMON/ETAGRD/
     *      XB,DEX,YB,DEY,
     *      COSETA(JM),XETA(IMAX),YETA(JM),
     *      ITMP(NPOINT), JTMP(NPOINT), KTMP(NPOINT)
     
      COMMON/ETAREA/CAREA(IMJM)
      
      COMMON/ETABND/IBND1,IBND2,JBND1,JBND2         
         
      COMMON/INNERP/IPN,INP(IMJM),IINP(IMJM),JINP(IMJM)      

      COMMON/TAB25A/IOA(25),WET(26),IOAC(16)
      
      COMMON/INNER2/IPN2,IPN3,INP2(IMJM),IIN2(IMJM),JIN2(IMJM),
     *                        INP3(IMJM)
      COMMON/CITERB/EPS,MAXIT          
                  
C-----------------------------------------------------------------------
                             C O M M O N /LOOPS/
     & KHLA  (JAM),KHHA  (JAM),KVLA  (JAM),KVHA  (JAM)
C
     &,KHL2  (JM),KHH2  (JM),KVL2  (JM),KVH2  (JM)
C
     &,LMH   (IMJM),LMV   (IMJM)   
C-----------------------------------------------------------------------
                             C O M M O N /MASKS/
     & HBM2  (IMJM),VBM2  (IMJM),VBM3  (IMJM),SM    (IMJM),SICE (IMJM)
C
     &,HTM   (IMJM,LM),VTM   (IMJM,LM)
C-----------------------------------------------------------------------
                             C O M M O N /VRBLS/
     & PD(IMJM),RES(IMJM),FIS(IMJM)
C
     &,U(IMJM,LM),V(IMJM,LM),T(IMJM,LM),Q(IMJM,LM)
C-----------------------------------------------------------------------
                              C O M M O N /MAPOT/
     & LSL,IXM,IYM
     &,ISHDE(49)
     &,R,PT,TSPH
     &,WBD,SBD,TLM0D,TPH0D,DLMD,DPHD,CMLD,DP30,X1P,Y1P
     &,DISLP,Z0SLP,SPL(LSM),ALSL(LSM)
     &,DETA(LM),AETA(LM),ETA(LP1),TSHDE(49)
     &,ERLAM0,CPHI0,SPHI0
     &,KUNIT,LUNIT,MUNIT,IUNIT,NUNIT
C------------------------------------------------------------------------
      COMMON/TM3COM/SURP8(8),STAX(10,20)       
      LOGICAL LFIL, UNDIS
    
      DATA PAI /3.14159265/
      DATA G /9.81/
      DATA RGAS /287.04/
      DATA ERAD /6370000./
      DATA CON /-999.0/
C      
C         S T A R T
      WRITE(66,*)'START TRPORT'
C/      WRITE(66,*)'IPN=',IPN
C/      WRITE(0,*)'START TRPORT'
C/      WRITE(0,*)'IPN=',IPN            
C      
      DO I=1,20
      STAX(1,I)=0.
      ENDDO
      
      PID180=PAI/180.
      CONV=PID180
C
C     XB, DEX
C     YB, DEY
C     CON        = -999.0
C
CX      CALL SETID3(XB,DEX,YB,DEY,IM,JM,CON)   ! See TRNETA                                             
C
C      DENSITY CALCULATION
C      &  ALFA CALCULATION
C
       DO K=1,KMAX
              
       DO I=1,IMJM
       
       IF(HTM(I,K+KMIN).EQ.1.)THEN
       
CX       PETA=PT+AETA(K+KMIN)*PD(I)*RES(I)   ! #############
CX       C1=PETA/RGAS 
             
       ALFA(I,K,1)=1.
       ALFA(I,K,2)=1.
       ELSE
       ALFA(I,K,1)=0.         ! CON                  !   -999.0
       ALFA(I,K,2)=0.         ! CON                  !   -999.0       
       ENDIF
       
       ENDDO
       ENDDO
             
      ILENG=IBND2-IBND1+1
      JLENG=JBND2-JBND1+1

      ICL2=0                        ! Debuging control      
      
      CX=1./(ERAD*PID180)
      
C/      WRITE(66,*)'IPN2=',IPN2,' IPN3=',IPN3
C/      WRITE(0,*)'IPN2=',IPN2,' IPN3=',IPN3
                
      IXYZ=0
CDB      IDEBUG=0
                  
                  
C/      WRITE(0,*)'START JP and KL Loops'
                  
      DO1000 JP=1,IPN
C/      WRITE(0,*)'JP=',JP      
      
CDB      NUMPOI=JP                !DB
      
      IH=INP(JP)      
      
      IVE=IH
      IVW=IVE-1
      IVN=IH+IM-1  
      IVS=IH-IM    

      I=IINP(JP)
      J=JINP(JP)
       
      CY=CX/COSETA(J)      
      
      XA=XB+(I-1)*DEX
      YA=YB+(J-1)*DEY         
      TOPA=TOPG(IH) 
      
CX      CALL SAMPLX(FLOAT(IH),STAX(1,9))
      
      DO3000 KL=1,KMAX
      
C/      WRITE(0,*)'KL=',KL
               
      KIT=0
      NSURP8=0
            
      IF(HTM(IH,KL+KMIN).EQ.0.)THEN 
C --------------------------------
C     UNDER SURFACE
C --------------------------------
C/      WRITE(0,*)'UNDER SURFACE'
      
      QQ(IH,K,2)=0.     
      ELSE
C --------------------------------
C     OVER SURFACE
C --------------------------------
C/      WRITE(0,*)'OVER SURFACE'
            
      QA=QQ(IH,KL,1) 
             
      ZA=AETA(KL+KMIN)
      WAR=WW(IH,KL,2)*HTM(IH,KL+KMIN)
      
      UA=VTM(IVE,KL+KMIN)*UU(IVE,KL,2)+
     *   VTM(IVW,KL+KMIN)*UU(IVW,KL,2)+
     *   VTM(IVS,KL+KMIN)*UU(IVS,KL,2)+
     *   VTM(IVN,KL+KMIN)*UU(IVN,KL,2)       
      
      VA=VTM(IVE,KL+KMIN)*VV(IVE,KL,2)+
     *   VTM(IVW,KL+KMIN)*VV(IVW,KL,2)+
     *   VTM(IVS,KL+KMIN)*VV(IVS,KL,2)+
     *   VTM(IVN,KL+KMIN)*VV(IVN,KL,2)
     
      IW=INT(VTM(IVE,KL+KMIN)+VTM(IVW,KL+KMIN)+
     *       VTM(IVN,KL+KMIN)+VTM(IVS,KL+KMIN))+1
      
      UA=UA*WET(IW)
      VA=VA*WET(IW)
      
      X1=XA-CX*UA*DT
      Y1=YA-CY*VA*DT
      Z1=ZA-CW*WAR*DT
      
      QQH(IH,KL)=0.              !  21/12/98
      QQL(IH,KL)=0.              !  21/12/98

C/      WRITE(0,*)'START SAMPLX'
      
      ONE=1.
CX      CALL SAMPLX(ONE,STAX(1,1))
CX      CALL SAMPLX(QA,STAX(1,2))
CX      CALL SAMPLX(UA,STAX(1,3))
CX      CALL SAMPLX(VA,STAX(1,4))
C/      WRITE(0,*)'IH =',IH,' KL=',KL       
C/      WRITE(0,*)'START LATERAL BOUNDARY FILTR'
                 
C      
C         LATERAL BOUNDARY FILTR
C
      CALL BOUNDT(X1,Y1,XETA(IBND1),ILENG,
     * YETA(JBND1),JLENG,IBOUND)
C/      WRITE(0,*)'IH =',IH,' KL=',KL      
C/       WRITE(0,*)'IBOUND=',IBOUND
      
      IF(IBOUND.NE.0)THEN      
      Q2=0.     
      QHP=Q2
      QLP=Q2
      QHP2=Q2
      NSURP8=0            
      GOTO999        
      ENDIF
C/      WRITE(0,*)'IH =',IH,' KL=',KL      
C/      WRITE(0,*)'START Z  FILTER'
C
C                Z  FILTER
C
      IF(Z1.LT.AETA(KTOP))THEN      !   Change 7/02/99
      Q2=0.                         !           /
      QHP=Q2                        !          /
      QLP=Q2                        !         /
      QHP2=Q2
      NSURP8=0                      !        / 
      GOTO999                       !       / 
      ENDIF                         !   Change 7/02/99  
                  
      Z1=AMAX1(Z1,AETA(KTOP))
      Z1=AMIN1(Z1,AETA(LM))
                 
C     ---------------------------------
C             X,Y, and Z  FILTER
C     ---------------------------------
      CALL LOCATE(AETA(KTOP),KMAX,Z1,LEVTOP) 
      
      IF(LEVTOP.LT.1)LEVTOP=1
      IF(LEVTOP.GE.KMAX)LEVTOP=KMAX-1  

      LEVBOT=LEVTOP+1        

C     ================= 
      
      KCELL=INDEXC(X1,Y1) 

C/      WRITE(0,*)'KCELL=',KCELL

      IF(LEVBOT.GT.LMH(KCELL))THEN
C               UNDER

C/      WRITE(0,*)'UNDER'
      NSURP8=0
      Q2=0.     
      QHP=Q2
      QLP=Q2
      QHP2=Q2
      GOTO999 
      ELSE
C               OVER
      ENDIF     
C     ------------------------------------      
C
C     FILTRATION OF SMALL VALUES (Q.LT. QMIN)
C
C     ------------------------------------
C/      WRITE(0,*)'START FILTRATION OF SMALL ...'
      
      CALL CFILT2
     * (X1,Y1,LEVTOP,LEVBOT,QQ(1,1,1),
     * HTM(1,KTOP),KMAX,QMIN,UNDIS)

      LFIL=.FALSE.

C/      WRITE(0,*)'UNDIS=',UNDIS

      IF(UNDIS)THEN
      Q2=QQ(IH,KL,1)
      Q22=QQ(IH,KL,1)
      DIF=0.
      QHP=Q2
      QLP=Q2
      QHP2=Q2
      NSURP8=0           
      GOTO999
      ENDIF      

      X0=X1
      Y0=Y1
      Z0=Z1
      X2=X1
      Y2=Y1
      Z2=Z1
      KIT=1      
      
      DQA=0. 
      
C/      WRITE(0,*)'START ITERATIONS'
                   
C     ----------------------------  ITERATIONS
C     ITERATIONS
C     ----------------------------  ITERATIONS
10    CONTINUE

      CALL BOUNDT(X2,Y2,XETA(IBND1),ILENG,
     * YETA(JBND1),JLENG,IBOUND)
      
      IF(IBOUND.NE.0)THEN      
      Q2=0.     
      QHP=Q2
      QLP=Q2
      QHP2=Q2
      NSURP8=0                  
      GOTO999        
      ENDIF
            
      Z2=AMAX1(Z2,AETA(KTOP))
      Z2=AMIN1(Z2,AETA(LM))
                             
C     ---------------------------------
C             X,Y, and Z  FILTER
C     ---------------------------------
      CALL LOCATE(AETA(KTOP),KMAX,Z2,LEVTOP) 
      
      IF(LEVTOP.LT.1)LEVTOP=1
      IF(LEVTOP.GE.KMAX)LEVTOP=KMAX-1  

      LEVBOT=LEVTOP+1        

C     ================= 
      
      KCELL=INDEXC(X2,Y2) 
       
      IF(LEVBOT.GT.LMH(KCELL))THEN
C               UNDER
      Q2=0.     
      QHP=Q2
      QLP=Q2
      QHP2=Q2
      NSURP8=0
      GOTO999 
      ELSE
C               OVER
      ENDIF       
C
C     INTERPOLATION FOR (X2,Y2,Z2) POINT
C 

C/      WRITE(0,*)
C/     *'START INTERPOLATION FOR (X2,Y2,Z2)'
 
C ............... UU(1,1,1)      
       CALL VUBI3D(X2,Y2,Z2,AETA(KTOP),KMAX,  
     * UU(1,1,1),VTM(1,KTOP),U2,HU2)
C ............... VV(1,1,1)
       CALL VUBI3D(X2,Y2,Z2,AETA(KTOP),KMAX,  
     * VV(1,1,1),VTM(1,KTOP),V2,HV2)
C ............... WW(1,1,1)
       CALL CUBI3D(X2,Y2,Z2,AETA(KTOP),KMAX,  
     * WW(1,1,1),HTM(1,KTOP),W2,HW2)

      IF(HU2*HV2*HW2.EQ.1)THEN   
      X2=XA-.5*CX*(UA+U2)*DT
      Y2=YA-.5*CY*(VA+V2)*DT
      Z2=ZA-.5*CW*(WAR+W2)*DT
                
      ELSE
      X2=X1
      Y2=Y1
      Z2=Z1    
      ENDIF

      D=DIST(X1,Y1,X2,Y2)
      IF(D.GT.EPS)THEN
      X1=X2
      Y1=Y2
      KIT=KIT+1
           IF(KIT.GE.MAXIT)THEN
           X2=X0
           Y2=Y0
           ELSE
           GOTO10
           ENDIF
      ENDIF
      DA2=DIST(XA,YA,X2,Y2)
C     ----------------------------  END OF ITERATIONS
11    CONTINUE
C     -----------------------------------
C
C                Z  FILTER
C      
      IF(Z2.GT.AETA(LM))THEN
      Z2=AETA(LM)
      ELSE
      IF(Z2.LT.AETA(KTOP))THEN
      Z2=AETA(KTOP)
      ENDIF
      ENDIF 
C     -----------------------------------     

      CALL CUBI3D(X2,Y2,Z2,AETA(KTOP),KMAX,  
     * QQ(1,1,1),HTM(1,KTOP),QHP2,HHP2)
      
      QHP=QHP2*HHP2
       
CX      CALL SAMPLX(HHP2,STAX(1,5))
CX      CALL SAMPLX(QHP,STAX(1,6))
            
      IXYZ=0 
     
      CALL LI3DSU(X2,Y2,Z2,AETA(KTOP),KMAX,
     * QQ(1,1,1),HTM(1,KTOP),SURP8,NSURP8,QLP,HQLP)
     
CX      CALL SAMPLX(QLP,STAX(1,7)) 
CX      CALL SAMPLX(HQLP,STAX(1,8))           
           
      IF(HQLP.NE.1.)QLP=0.
C ----------------------------- 
C     QUASI-MONOTONIC ALGORITHM
C -----------------------------
      IF(NSURP8.NE.0)THEN
      QMN=SURP8(1)
      QMX=SURP8(1)
      DO I=1,NSURP8
      QMN=AMIN1(QMN,SURP8(I))
      QMX=AMAX1(QMX,SURP8(I))
      ENDDO
      
      QHP=AMAX1(QHP,QMN)
      QHP=AMIN1(QHP,QMX)
      ELSE
      QHP=0.
      ENDIF
CX            O.K.      
CX      IF(QHP.LT.0.)THEN
CX      QHP=0.
CX      ENDIF
C ----------------------------- 
               
      QQH(IH,KL)=QHP      
      QQL(IH,KL)=QLP 
    
C     --------------
C     Internal point
C     --------------

C     ---------------------------------
C            BOUNDARY?
C     ---------------------------------
 999    CONTINUE
 
C/      WRITE(0,*)'999'
C/      WRITE(0,*)'NSURP8=',NSURP8      
C     ---------------------------------sem diffusion  
      Q3=QHP
      Q23=QLP
C     --------------------------------- 
      IF(QHP.LT.0..AND.QHP.NE.CON)THEN
      QHP=0.
      ENDIF
      
C/      WRITE(0,*)'QHP=',QHP 
C/      WRITE(0,*)'IH =',IH,' KL=',KL
        
      QQ(IH,KL,2)=QHP    

C     ---------------------------------------------------------
C                    Here I can calculate ALFA
C     ---------------------------------------------------------
      
C/      WRITE(0,*)'NSURP8=',NSURP8
      
          IF(NSURP8.NE.0)THEN
C/          WRITE(0,*)'SURP8=',SURP8
C
      CALL ALFMAX(QHP,QLP,QMP,SURP8,NSURP8,ALFA(IH,KL,2))
C
C/          WRITE(0,*)'QHP=',QHP
C/          WRITE(0,*)'QLP=',QLP 
C/          WRITE(0,*)'QMP=',QMP
          ENDIF     
C     ---------------------------------------------------------
      ENDIF    ! ONE POINT CALCULATION OVER SURFACE END

CX      CALL SAMPLX(QHP,STAX(1,8))
CX      CALL SAMPLX(QLP,STAX(1,9))
      
C/      WRITE(0,*)'End of loop KL'
                    
 3000     CONTINUE               !   LEVEL LOOP
 1000     CONTINUE           !   J LOOP
 
C/      CALL WRISAM(0,'N p. ',STAX(1,1))
C/      CALL WRISAM(0,'QHP  ',STAX(1,6))
            
CX      CALL WRISAM(66,'N p. ',STAX(1,1))
CX      CALL WRISAM(66,'QA   ',STAX(1,2)) 
CX      CALL WRISAM(66,'UA   ',STAX(1,3)) 
CX      CALL WRISAM(66,'VA   ',STAX(1,4))
CX      CALL WRISAM(66,'HHP2 ',STAX(1,5))
CX      CALL WRISAM(66,'QHP  ',STAX(1,6))
CX      CALL WRISAM(66,'QLP  ',STAX(1,7))
CX      CALL WRISAM(66,'2.QHP',STAX(1,8))
CX      CALL WRISAM(66,'2.QLP',STAX(1,9))            
      
      
      CALL EXTRAC(QQ(1,1,1),IMJM*KMAX,AMEAN,AMIN,AMAX,0.,NULL)
      WRITE(66,*)'Q(K  ):  AMEAN,AMIN,AMAX,NULL=',
     * AMEAN,AMIN,AMAX,NULL
      CALL EXTRAC(QQ(1,1,2),IMJM*KMAX,AMEAN,AMIN,AMAX,0.,NULL)
      WRITE(66,*)'Q(K+1):  AMEAN,AMIN,AMAX,NULL=',
     * AMEAN,AMIN,AMAX,NULL 
      CALL EXTRAC(HTM,IMJM*KMAX,AMEAN,AMIN,AMAX,0.,NULL)
      WRITE(66,*)'HTM (0) :  AMEAN,AMIN,AMAX,NULL=',
     * AMEAN,AMIN,AMAX,NULL 
      CALL EXTRAC(HTM,IMJM*KMAX,AMEAN,AMIN,AMAX,1.,NULL)
      WRITE(66,*)'HTM (1) :  AMEAN,AMIN,AMAX,NULL=',
     * AMEAN,AMIN,AMAX,NULL                             
                               
      RETURN
   
      END
C& BOUNDT
      SUBROUTINE BOUNDT(XP,YP,X,M,Y,N,IB)
      DIMENSION X(M),Y(N)
      IB=0
      IF(XP.LT.X(1))IB=1000
      IF(XP.GT.X(M))IB=IB+100
      IF(YP.LT.Y(1))IB=IB+10
      IF(YP.GT.Y(N))IB=IB+1
      RETURN
      END
C& VFLUXS
      SUBROUTINE VFLUXS(NU,ISTEP,QFLUX,EFIREA)    !  ,BUFI,IRECF)
C
C     SURFACE FLUX
C
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)      
      PARAMETER (IMJM = IM*JM - JM/2)
      PARAMETER (LP1=LM+1)
      PARAMETER (JAM=6+2*(JM-10))      
      PARAMETER (IMAX=2*IM-1  , JMAX=JM)       
      PARAMETER (KKMAX=LM)
      PARAMETER (KMAX=KKMAX)
      PARAMETER (KMIN = LM-KKMAX)
      
                 
      DIMENSION QFLUX(IMJM,2),EFIREA(IMJM,4)   !  ,BUFI(IMJM,2)
      DIMENSION VFLX(14)
C-----------------------------------------------------------------------      
      COMMON/ETAVEG/IHVEG(IMJM)
C-----------------------------------------------------------------------
                             C O M M O N /LOOPS/
     & KHLA  (JAM),KHHA  (JAM),KVLA  (JAM),KVHA  (JAM)
C
     &,KHL2  (JM),KHH2  (JM),KVL2  (JM),KVH2  (JM)
C
     &,LMH   (IMJM),LMV   (IMJM)     
C------------------------------------------------------------------------
                              C O M M O N /MAPOT/
     & LSL,IXM,IYM
     &,ISHDE(49)
     &,R,PT,TSPH
     &,WBD,SBD,TLM0D,TPH0D,DLMD,DPHD,CMLD,DP30,X1P,Y1P
     &,DISLP,Z0SLP,SPL(LSM),ALSL(LSM)
     &,DETA(LM),AETA(LM),ETA(LP1),TSHDE(49)
     &,ERLAM0,CPHI0,SPHI0
     &,KUNIT,LUNIT,MUNIT,IUNIT,NUNIT
C-----------------------------------------------------------------------
                             C O M M O N /VRBLS/
     & PD(IMJM),RES(IMJM),FIS(IMJM)
C
     &,U(IMJM,LM),V(IMJM,LM),T(IMJM,LM),Q(IMJM,LM)
C-----------------------------------------------------------------------
      COMMON/ETATRN/
     * QQL(IMJM,KKMAX),QQH(IMJM,KKMAX),      
     * PS(IMJM,2), QQ(IMJM,KKMAX,2),  
     * UU(IMJM,KKMAX,2), VV(IMJM,KKMAX,2),
     * WW(IMJM,KKMAX,2), RO(IMJM,KKMAX,2),     
     * TT(IMJM,KKMAX,2), ZZ(IMJM,KKMAX+1,2),
     * ALFA(IMJM,KKMAX,2),
     * DIFQ(IMJM,KKMAX),
     * AKH(IMJM,KKMAX,2),AKHS(IMJM,2),
     * TOPG(IMJM)
     *,PREC(IMJM,2),ACPR(IMJM,2)
C-----------------------------------------------------------------------
C      
      COMMON/ETAREA/CAREA(IMJM)
C                 0  1 2 3 4 5 6 7 8 9 10 11 12 13 
      DATA VFLX/ 0.,                 ! 0 ocean
     *           0.206E-3,           ! 1 floresta
     *           0.580E-4,           ! 2
     *           0.580E-4,           ! 3     
     *           0.580E-4,           ! 4     
     *           0.580E-4,           ! 5       
     *           0.903E-6,           ! 6 cerrado 
     *           0.580E-4,           ! 7                
     *           0.580E-4,           ! 8      
     *           0.580E-4,           ! 9                
     *           0.580E-4,           ! 10
     *           0.580E-4,           ! 11      
     *           0.580E-4,           ! 12
     *           0.580E-4/           ! 13        
C           
      DATA RGAS/287.04/    
C-----------------------------------------------------------------------          
      IHR=MOD(ISTEP,30)
C
C                 QFLUX(*,1)
C      
C
C     EFIREA(J) in km**2
C     CAREA(J)  in m**2
C
C      C1=QFLX*1.E6  ! For converting EFIREA(J) to m**2
C
      C2=1.E12    ! For converting to ppm  and EFIREA(J) to m**2         
      DO J=1,IMJM
      IVT=IHVEG(J)
      
      IF(IVT.GT.13)THEN
      WRITE(66,*)'IVT=',IVT,' for ',J,' point'
      STOP      
      ENDIF
      
      QFLXV=VFLX(IVT+1)   !*0.5    !  equiv fire area change by 0.5 times
      QFLUX(J,1)=QFLXV*C2*EFIREA(J,1)/CAREA(J)
      EFIREA(J,2)=0.
      EFIREA(J,3)=0.
      EFIREA(J,4)=0.
      ENDDO      
C
C                 QFLUX(*,2)
C
       DO10 J=1,IMJM
       
       C=PD(J)*RES(J)
       K=LMH(J)
       IF(K.EQ.0)THEN
       QFLUX(J,2)=0.
       QFLUX(J,1)=0.
       ELSE
       PSUR=PT+AETA(K+KMIN)*C
       TEM=TT(J,K,2)
       ROP=PSUR/(RGAS*TEM)
       DZ=ZZ(J,K,2)-ZZ(J,K+1,2)
       QFLUX(J,2)=QFLUX(J,1)/(ROP*DZ)
       ENDIF
       
 10    CONTINUE
C      
      RETURN
      END       
C& SFLUXS
      SUBROUTINE SFLUXS(NU,ISTEP,QFLUX,QFLX,EFIREA,BUFI,IRECF)
C
C     SURFACE FLUX
C
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)      
      PARAMETER (IMJM = IM*JM - JM/2)
      PARAMETER (LP1=LM+1)
      PARAMETER (JAM=6+2*(JM-10))      
      PARAMETER (IMAX=2*IM-1  , JMAX=JM)       
      PARAMETER (KKMAX=LM)
      PARAMETER (KMAX=KKMAX)
      PARAMETER (KMIN = LM-KKMAX)
      
                 
      DIMENSION QFLUX(IMJM,2),EFIREA(IMJM,4),BUFI(IMJM,2)
C-----------------------------------------------------------------------
                             C O M M O N /LOOPS/
     & KHLA  (JAM),KHHA  (JAM),KVLA  (JAM),KVHA  (JAM)
C
     &,KHL2  (JM),KHH2  (JM),KVL2  (JM),KVH2  (JM)
C
     &,LMH   (IMJM),LMV   (IMJM)     
C------------------------------------------------------------------------
                              C O M M O N /MAPOT/
     & LSL,IXM,IYM
     &,ISHDE(49)
     &,R,PT,TSPH
     &,WBD,SBD,TLM0D,TPH0D,DLMD,DPHD,CMLD,DP30,X1P,Y1P
     &,DISLP,Z0SLP,SPL(LSM),ALSL(LSM)
     &,DETA(LM),AETA(LM),ETA(LP1),TSHDE(49)
     &,ERLAM0,CPHI0,SPHI0
     &,KUNIT,LUNIT,MUNIT,IUNIT,NUNIT
C-----------------------------------------------------------------------
                             C O M M O N /VRBLS/
     & PD(IMJM),RES(IMJM),FIS(IMJM)
C
     &,U(IMJM,LM),V(IMJM,LM),T(IMJM,LM),Q(IMJM,LM)
C-----------------------------------------------------------------------
      COMMON/ETATRN/
     * QQL(IMJM,KKMAX),QQH(IMJM,KKMAX),      
     * PS(IMJM,2), QQ(IMJM,KKMAX,2),  
     * UU(IMJM,KKMAX,2), VV(IMJM,KKMAX,2),
     * WW(IMJM,KKMAX,2), RO(IMJM,KKMAX,2),     
     * TT(IMJM,KKMAX,2), ZZ(IMJM,KKMAX+1,2),
     * ALFA(IMJM,KKMAX,2),
     * DIFQ(IMJM,KKMAX),
     * AKH(IMJM,KKMAX,2),AKHS(IMJM,2),
     * TOPG(IMJM)
     *,PREC(IMJM,2),ACPR(IMJM,2)
C-----------------------------------------------------------------------
C      
      COMMON/ETAREA/CAREA(IMJM)
C
      DATA RGAS/287.04/    
C-----------------------------------------------------------------------          
      IHR=MOD(ISTEP,30)
C        
C      CALL RDBUF(NU,ISTEP,IRECF,BUFI,EFIREA,IMJM) 
C          RDBUF(NU,ISTEP,IREC,BUFI,FIRE,IMJM)
C
C                 QFLUX(*,1)
C      
C
C     EFIREA(J) in km**2
C     CAREA(J)  in m**2
C
      C1=QFLX*1.E6  ! For converting EFIREA(J) to m**2
      C2=C1*1.E6    ! For converting to ppm          
      DO J=1,IMJM
      QFLUX(J,1)=C2*EFIREA(J,1)/CAREA(J)
      EFIREA(J,2)=0.
      EFIREA(J,3)=0.
      EFIREA(J,4)=0.
      ENDDO      
C
C                 QFLUX(*,2)
C
       DO10 J=1,IMJM
       
       C=PD(J)*RES(J)
       K=LMH(J)
       IF(K.EQ.0)THEN
       QFLUX(J,2)=0.
       QFLUX(J,1)=0.
       ELSE
       PSUR=PT+AETA(K+KMIN)*C
       TEM=TT(J,K,2)
       ROP=PSUR/(RGAS*TEM)
       DZ=ZZ(J,K,2)-ZZ(J,K+1,2)
       QFLUX(J,2)=QFLUX(J,1)/(ROP*DZ)
       ENDIF
       
 10    CONTINUE
C      
      RETURN
      END
C& RDBUF
      SUBROUTINE RDBUF(NU,ISTEP,IREC,BUFI,FIRE,IMJM)
      DIMENSION BUFI(IMJM,2),FIRE(IMJM)
      DIMENSION INX(30), ALPHA(4)
      DATA INX/1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,
     *         3,3,3,3,3,3,3,4,4,4,4,4,4,4,4/
      DATA ALPHA/0.47,0.53,0.56,0.44/
         
      IDAY   = (ISTEP-1)/30+1      !             day
      TIME   = (IDAY-1)*24.*60.    !      min    time
      
      ISTED  = MOD(ISTEP-1,30)+1   !             inner day step
      TIMED  = ISTED*48.           !      min    inner time
             
          IF(ISTED.EQ.1.OR.ISTED.EQ.16)THEN       !  new
C                                   15            !  old   	 
	 READ(NU,rec=IREC)(BUFI(J,1),J=1,IMJM)       
	 	            
      WRITE(66,*)
     *'FIRE AREA OF IDAY=',IDAY,' HOURE=',TIME,' IREC=',IREC,' /#' 
	 IREC=IREC+1
	 
	 IF(IREC.EQ.3)IREC=1          !  new   Next day is percisten (second)
	                      	 	 
	 ENDIF
	 
             IA=INX(ISTED)              !  new
             ALP=ALPHA(IA)              !  new
             	 
             DO J=1,IMJM
             FIRE(J)=ALP*BUFI(J,1)      !  new            
             ENDDO		  	 	                        	  
      RETURN   
      END       
C& WRSUR
      SUBROUTINE WRSUR(NU,IREC,AA,HZM,BB,CON)           
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)
      PARAMETER (LP1=LM+1)            
      PARAMETER (IMJM = IM*JM - JM/2)
      PARAMETER (JM2=JM/2+1,IMT=IM*2-1)      
      PARAMETER (NPOINT=25)
      PARAMETER (IMAX=51, JMAX=51)      
C
      DIMENSION AA(IMJM),HZM(IMJM),BB(IMAX,JMAX)
      COMMON/ETAGRD/
     * XB,DEX,YB,DEY,
     * COSETA(JM),XETA(IMT),YETA(JM),
     * ITMP(NPOINT),JTMP(NPOINT),KTMP(NPOINT)
C----------------------------------------------
       C O M M O N /MAPOT/
     & LSL,IXM,IYM
     &,ISHDE(49)
     &,R,PT,TSPH
     &,WBD,SBD,TLM0D,TPH0D,DLMD,
     & DPHD,CMLD,DP30,X1P,Y1P
     &,DISLP,Z0SLP,SPL(LSM),ALSL(LSM)
     &,DETA(LM),AETA(LM),ETA(LP1),TSHDE(49)
     &,ERLAM0,CPHI0,SPHI0
     &,KUNIT,LUNIT,MUNIT,IUNIT,NUNIT
C-----------------------------------------------
       C O M M O N /VRBLS/
     & PD(IMJM),RES(IMJM),FIS(IMJM)
C
     &,U(IMJM,LM),V(IMJM,LM),
     & T(IMJM,LM),Q(IMJM,LM)
C---------------------------------------------          
      DATA ERAD /6370000./
C
      ICL=0
C      
       RADDEG=3.14159265/180.
       DEG=180./3.14159265
C
       DPHI=DPHD*RADDEG   !   in radian 
       DLAM=DPHD*RADDEG   !   in radian 
       
       GLON0=TLM0D*RADDEG
       GLAT0=TPH0D*RADDEG
       
       DLON=1.0
       DLAT=1.0
       
       DO J=1,JMAX
       YJ=-40.0+(J-1)*DLAT
       GLAT=YJ*RADDEG           
       DO I=1,IMAX
       XI=-80.0+(I-1)*DLON
       GLON=XI*RADDEG 
        
       CALL ToETA(GLON,GLAT,
     * ELON,ELAT,GLON0,GLAT0)
      
       XE=ELON*DEG
       YE=ELAT*DEG 

       IF(
     * XE.GE.XETA(3).AND.XE.LE.XETA(IMT-2)
     *                .AND.
     * YE.GE.YETA(3).AND.YE.LE.YETA(JM-2)
     *   )THEN
       
       CALL AV25P(XE,YE,AA,HZM,IMJM,CON,QP)     
            
       BB(I,J)=QP
                       
       ELSE
C
C      OUT MODEL AREA
C       
       BB(I,J)=CON 
               
       ENDIF           


       ENDDO     ! I loop       
       ENDDO     ! J loop 
C       
C     FILTR
C
      DO J=1,JMAX
      DO I=1,IMAX
      IF(ABS(BB(I,J)).LT.1.E-20)
     * BB(I,J)=0.
      ENDDO
      ENDDO
         
      IREC=IREC+1
      WRITE(NU,rec=IREC)
     * ((BB(I,J),I=1,IMAX),J=1,JMAX)
 
       RETURN
       END
C& SUBTRA
      SUBROUTINE SUBTRA(A,B,M,C)
C     CALL SUBTRA(APRC(1,2),APRC(1,1),IMJM,DPREC)
      DIMENSION A(M), B(M), C(M)
      DO J=1,M
      C(J)=A(J)-B(J)
      ENDDO
      RETURN
      END
C& WRH6VR          !!!!!!!!!!!!!!!!!! nvp
      SUBROUTINE WRH6VR(NU,IREC,QQ,RO,ZZ,AU,AV,
     * HTM,HZM,VTM,LMH,PLEV1,PLEV2,BB,BU,BV,WA,KMAX,CON)
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)
      PARAMETER (LP1=LM+1)            
      PARAMETER (IMJM = IM*JM - JM/2)
      PARAMETER (JM2=JM/2+1,IMT=IM*2-1)      
      PARAMETER (IMAX=51, JMAX=51)
      DIMENSION HTM(IMJM,KMAX),HZM(IMJM,KMAX+1),
     * VTM(IMJM,KMAX),LMH(IMJM)      
      DIMENSION QQ(IMJM,KMAX),RO(IMJM,KMAX),ZZ(IMJM,KMAX+1),
     * AU(IMJM,KMAX),AV(IMJM,KMAX),
     * BB(IMAX,JMAX,6),
     * BU(IMAX,JMAX,6),BV(IMAX,JMAX,6),WA(IMJM,KMAX),
     * PLEV1(KMAX),PLEV2(KMAX+1)
      COMMON/ETAREA/CAREA(IMJM)
C ----------------------------------------------
CX      PARAMETER (LEVS = 6)
CX      PARAMETER (NEL3 = 4)
CX      PARAMETER (NEL2 = 2)
CX      PARAMETER (NSTEP = 9)
C -----------------------------------------------
C     0   12  24  36  48  60  72  84  96  hours
C     x---o---x---o---x---o---x---o---x
C     1   2   3   4   5   6   7   8   9   step
C -----------------------------------------------
CX      PARAMETER (MAXREC=(LEVS*NEL3+NEL2)*NSTEP)      
C ----------------------------------------------
CX      IREST=MOD(IREC,MAXREC)
CX      IF(IREST.EQ.0)THEN
CX      
CX      WRITE(66,*)'FIRST HISTORY: IREC=',IREC,' MAXREC=',MAXREC,' #/'
CX      
CX      IREC=0      
CX      ELSE
CX      IREC=IREST        
CX      ENDIF                  

      WRITE(66,*)'START HISTORY RECORD=',IREC,' #/'
      WRITE(0 ,*)'START HISTORY RECORD=',IREC,' #/'
CX         LOOKHP(NU,NAME,A,M,KMAX,HTM)
CX      CALL LOOKHP(66,'RO',RO,IMJM,KMAX,HTM)
CX      CALL LOOKHP(66,'QQ',QQ,IMJM,KMAX,HTM)
       
      DO K=1,KMAX
      DO J=1,IMJM
      
CX      IF(QQ(J,K).NE.CON.AND.RO(J,K).NE.CON)THEN
CX      WA(J,K)=QQ(J,K)*RO(J,K)*1000.0      !  mikro g / m*3
CX      ELSE
CX      WA(J,K)=CON
CX      ENDI
c      C1=1./(CAREA(J)*(ZZ(J,K)-ZZ(J,K+1)))
c      WA(J,K)=QQ(J,K)*RO(J,K)*1000.0*HTM(J,K)*C1      !  mikro g / m*3
      WA(J,K)=QQ(J,K) 
      ENDDO
      ENDDO
      
      CALL LOOKHP(66,'rho_a',WA,IMJM,KMAX,HTM)           
                             
      CALL WROUT(NU,IREC,WA,HTM,HZM,
     * PLEV1,KMAX,BB,CON)
      CALL EXTRAC(BB(1,1,3),IMAX*JMAX,
     * AMEAN,AMIN,AMAX,-999.0,NCON)
      WRITE(0,*)'QQ(1,1,3)=',AMEAN,AMIN,AMAX,-999.0,NCON     
      WRITE(0,*)'FINISH WROUT(...,QQ... mikro g / m*3'
      
      DO J=1,IMJM
      KSURF=LMH(J) 
      QSURF=WA(J,KSURF)         
      WA(J,1)=QSURF
      ENDDO 
      
      CALL LOOKHP(66,'QS',WA,IMJM,1,HZM)            
      
      CALL WROUT(NU,IREC,ZZ,HZM,HZM,
     * PLEV2,KMAX+1,BB,CON)
      CALL EXTRAC(BB(1,1,3),IMAX*JMAX,
     * AMEAN,AMIN,AMAX,-999.0,NCON)
      WRITE(0,*)'ZZ(1,1,3)=',AMEAN,AMIN,AMAX,-999.0,NCON
            
      WRITE(0,*)'FINISH WROUT(...,ZZ'     
     
      CALL WROUV(NU,IREC,AU,AV,VTM,HZM,PLEV1,KMAX,
     *              BU,BV,CON)
     
      CALL WRSUR(NU,IREC,WA,HZM,BB,CON)
           
      CALL EXTRAC(BU(1,1,3),IMAX*JMAX,
     * AMEAN,AMIN,AMAX,-999.0,NCON)
      WRITE(0,*)'UU(1,1,3)=',AMEAN,AMIN,AMAX,-999.0,NCON 
      CALL EXTRAC(BV(1,1,3),IMAX*JMAX,
     * AMEAN,AMIN,AMAX,-999.0,NCON)
      WRITE(0,*)'VV(1,1,3)=',AMEAN,AMIN,AMAX,-999.0,NCON           
      WRITE(0,*)'FINISH WROUV(...,AU,AV'
      
      CALL EXTRAC(BB,IMAX*JMAX,
     * AMEAN,AMIN,AMAX,-999.0,NCON)
      WRITE(0,*)'QSURF=',AMEAN,AMIN,AMAX,-999.0,NCON           
      WRITE(0,*)'FINISH WRSUR(...,WA'      
           
      WRITE(66,*)'FINISH HISTORY RECORD=',IREC,' #/'
      WRITE( 0,*)'FINISH HISTORY RECORD=',IREC,' #/'

CX      STOP

      RETURN
      END         
C& LOOKHP 
      SUBROUTINE LOOKHP(NU,NAME,A,M,KMAX,HTM)
      DIMENSION A(M,KMAX),HTM(M,KMAX) 
      CHARACTER*(*) NAME
      WRITE(NU,*)'ARRAY ',NAME
      DO K=1,KMAX

      S1=0.
      N1=0     !   points with value
      N2=0     !   points without value
      DO J=1,M
        IF(HTM(J,K).EQ.1.)THEN
            IF(N1.EQ.0)THEN
            AMN=A(J,K)
            AMX=A(J,K)
            ENDIF
        S1=S1+A(J,K)            
        IF(A(J,K).LT.AMN)AMN=A(J,K)
        IF(A(J,K).GT.AMX)AMX=A(J,K)
        N1=N1+1
        ELSE
        N2=N2+1
        ENDIF
      ENDDO    !   Loop M      
      IF(N1.NE.0)THEN
      AMEAN=S1/FLOAT(N1)
      ELSE
      AMEAN=0.
      ENDIF
      AMIN=AMN
      AMAX=AMX      
      NCON=N2      
      WRITE(NU,100)K,AMEAN,AMIN,AMAX,NCON
      ENDDO    !   Loop KMAX      
      RETURN
100   FORMAT(I3,1X,3(1X,G12.4),2X,I4)
      END
C& SFLXD
      SUBROUTINE SFLXD(NU,NAME,QFLUX,EFIREA,CAREA,HZM,
     *   IMJM,QFLX,DTT)
      DIMENSION QFLUX(IMJM),EFIREA(IMJM),CAREA(IMJM),HZM(IMJM)
      CHARACTER*(*) NAME
C
      WRITE(NU,*)'  '
      WRITE(NU,*)'Distribution Source Diagnosis ',NAME
            
      CALL EXTRAC(QFLUX,IMJM,AMEAN,AMIN,AMAX,0.,NULL)      
      CONV1=1.E4*86400.*1.E-6   !  Convert from (m**2 s) to (ha day)

      WRITE(NU,*)'AMEAN=',AMEAN*1.E-6,
     *' AMIN=',AMIN*1.E-6,' AMAX=',AMAX*1.E-6,' kg/m**2/s'
      WRITE(NU,*)'NULL=',NULL,' #/'
      WRITE(NU,*)'AMEAN=',AMEAN*CONV1,
     *' AMIN=',AMIN*CONV1,' AMAX=',AMAX*CONV1,' kg/ha/day'
            
      CALL SCLINP(QFLUX,CAREA,HZM,IMJM,1.E-6,SUMFLX)
      WRITE(66,*)'1] DISTRIBUTED FLUX FROM SURFACE=',SUMFLX/1000.,
     * ' t/s'
     
      WRITE(NU,*)'2] DISTRIBUTED FLUX FROM SURFACE=',
     * DTT*SUMFLX/1000.0,' t/step  1 step=',DTT,' sec' 
         
      QFLUX0=DTT*QFLX*1.E+6
      CALL SCLINP(EFIREA,HZM,HZM,IMJM,QFLUX0,SUMFX1)
      WRITE(NU,*)'3] DISTRIBUTED FLUX FROM SURFACE=',
     * SUMFX1/1000.,' t/step'
     
      CALL SCLINP(EFIREA,HZM,HZM,IMJM,1.,SUMFX2)
      WRITE(NU,*)'4] FIRE AREA=',SUMFX2,' km**2 ',SUMFX2*1.E6,' m**2'
      
      RETURN
      END          
C& CLMDOM
      SUBROUTINE CLMDOM(Q,H,RES,PD,DETA,KMAX,
     *      COLCON,COLAIR)
C
C     COLUMN DENSITY OF MASS
C     COLCON   - column mass density of q in kg x 1.E+6
C     COLAIR   - column mass density of air in kg x 1.E+6
C
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)      
      PARAMETER (IMJM = IM*JM - JM/2)
      PARAMETER (LP1=LM+1)
      PARAMETER (IMAX=2*IM-1  , JMAX=JM)     

      DIMENSION Q(IMJM,KMAX),H(IMJM,KMAX),RES(IMJM),DETA(KMAX),
     *          PD(IMJM),
     *          COLCON(IMJM), COLAIR(IMJM)
      COMMON/INNERP/IPN,INP(IMJM),IINP(IMJM),JINP(IMJM) 
C
      DATA G /9.81/
C      
      DO J=1,IPN
      JP=INP(J)      

         AV=0.
         AM=0.

         C=RES(JP)*PD(JP)
               
         DO K=1,KMAX
         DM=H(JP,K)*DETA(K)*C
         AM=AM+DM  
         AV=AV+Q(JP,K)*DM
         ENDDO

      COLCON(JP)=1.E-6*AV/G   !  kg/m^2
      COLAIR(JP)=AM/G         !  kg/m^2  
      
      ENDDO

      RETURN
      END       
C& SITIOS
      SUBROUTINE SITIOS(NU,WA)   !   
      PARAMETER (NS=8)
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)
      PARAMETER (LP1=LM+1)            
      PARAMETER (IMJM = IM*JM - JM/2)
      PARAMETER (JM2=JM/2+1,IMT=IM*2-1)
      PARAMETER (IMAX=IMT)      
      PARAMETER (NPOINT=25)
      PARAMETER (KKMAX=LM)
      PARAMETER (KMAX=KKMAX)
      PARAMETER (KMIN = LM-KKMAX)
      PARAMETER (KTOP = KMIN + 1)                  
C
      DIMENSION WA(IMJM,KKMAX) 
C                                
      COMMON/ETAGRD/
     * XB,DEX,YB,DEY,
     * COSETA(JM),XETA(IMT),YETA(JM),
     * ITMP(NPOINT),JTMP(NPOINT),KTMP(NPOINT)
C----------------------------------------------
       C O M M O N /MAPOT/
     & LSL,IXM,IYM
     &,ISHDE(49)
     &,R,PT,TSPH
     &,WBD,SBD,TLM0D,TPH0D,DLMD,
     & DPHD,CMLD,DP30,X1P,Y1P
     &,DISLP,Z0SLP,SPL(LSM),ALSL(LSM)
     &,DETA(LM),AETA(LM),ETA(LP1),TSHDE(49)
     &,ERLAM0,CPHI0,SPHI0
     &,KUNIT,LUNIT,MUNIT,IUNIT,NUNIT
C-----------------------------------------------
                             C O M M O N /MASKS/
     & HBM2  (IMJM),VBM2  (IMJM),VBM3  (IMJM),SM    (IMJM),SICE (IMJM)
C
     &,HTM   (IMJM,LM),VTM   (IMJM,LM)  
C ---------------------------------------------------------------------
       C O M M O N /VRBLS/
     & PD(IMJM),RES(IMJM),FIS(IMJM)
C
     &,U(IMJM,LM),V(IMJM,LM),
     & T(IMJM,LM),Q(IMJM,LM)
C---------------------------------------------           
      COMMON/ETATRN/
     * QQL(IMJM,KKMAX),QQH(IMJM,KKMAX),      
     * PS(IMJM,2), QQ(IMJM,KKMAX,2),  
     * UU(IMJM,KKMAX,2), VV(IMJM,KKMAX,2),
     * WW(IMJM,KKMAX,2), RO(IMJM,KKMAX,2),     
     * TT(IMJM,KKMAX,2), ZZ(IMJM,KKMAX+1,2),
     * ALFA(IMJM,KKMAX,2),
     * DIFQ(IMJM,KKMAX),
     * AKH(IMJM,KKMAX,2),AKHS(IMJM,2),
     * TOPG(IMJM)
     *,PREC(IMJM,2),ACPR(IMJM,2)
C
      COMMON/CBT1T2/HINTG(IMJM), HINN(IMJM,KMAX),
     *              HZM(IMJM,KMAX+1),BUFI(IMJM,2),      
     *              QFLUX(IMJM,2) ,EFIREA(IMJM,4),
     *              QSOUR(10),ISOUR(10),JSOUR(10),
     *              XSOUR(10),YSOUR(10),
     *              QFIL(IMJM,KKMAX,2),     
     *              IDATCP(3)     
C ---------------------------------------------------------------------    
C |      
C |                         Latitude  Longitude  Altitude
C |                     (S)       (W)        (m)
C |                          
C | Fazenda Boa Sorte  5.17      48.75       170
C | Fazenda Dimona     2.32      60.32       120
C | Fazenda N.S.
C | Aparecida         10.75      62.37       220
C | Mina Potossi       9.27      62.87        80
C | Reserva Ducke      2.95      59.95        80
C | Reserva Jarú      10.08      61.92       120
C | Reserva Vale do
C | Rio Doce           5.75      49.17       150
C |
      DIMENSION X(NS),Y(NS),ALT(NS),WA1(NS),HMSK(NS),OD(NS),
     *          WC(NS),WPP(NS),WP(NS),WCD(LM,NS)
      CHARACTER*27 SITE(NS)
      DATA  X/-48.75,-60.32,-62.37,-62.87,-59.95,-61.92,-49.17,0./        
      DATA  Y/ -5.17, -2.32,-10.75, -9.27, -2.95,-10.08, -5.75,0./    
      DATA ALT/170.,  120.,  220.,   80.,   80.,  120.,  150., 0. /
      DATA SITE/     
     *'Fazenda Boa Sorte          ',     
     *'Fazenda Dimona Fazenda N.S.',
     *'Aparecida                  ',
     *'Mina Potossi               ',
     *'Reserva Ducke              ',
     *'Reserva Jarú               ',
     *'Reserva Vale do Rio Doce   ',
     *'Point with maximal concent.'/ 
C      123456789 123456789 1234567' 
C23456
      DATA ERAD /6370000./
      DATA G /9.81/      
C ------------------------------------
C  
C            CALCULATION
C
C ------------------------------------
      QMAX=0.          
      JMAX=1
      
      DO J=1,IMJM
      WA(J,2)=(PS(J,1)-PT)*RES(J)      !  N/m^2
      ENDDO
  
      DO J=1,IMJM
      PPS=(PS(J,1)-PT)*RES(J)
      C1=1.E-3*PPS/G
      S=0.
      
      DO K=1,KKMAX                          
      S=S+C1*QQ(J,K,1)*DETA(K)*HTM(J,K)
      ENDDO
           
      DO K=1,KKMAX        
          IF(QQ(J,K,1).GT.QMAX)THEN
          JMAX=J
          QMAX=QQ(J,K,1)
          ENDIF         
      ENDDO                              !  g/m^2
      
      WA(J,1)=S      
      
      ENDDO  !   END of PRESSURE WA(*,2) and COLUMN MASS WA(*,1)
C
      ICL=0
C      
       RADDEG=3.14159265/180.
       DEG=180./3.14159265
C
       DPHI=DPHD*RADDEG   !   in radian 
       DLAM=DPHD*RADDEG   !   in radian 
       
       GLON0=TLM0D*RADDEG
       GLAT0=TPH0D*RADDEG

       
      DO J=1,NS-1
               
      YJ=Y(J)
      GLAT=YJ*RADDEG           
      XI=X(J)
      GLON=XI*RADDEG 
            
      CALL ToETA(GLON,GLAT,
     *ELON,ELAT,GLON0,GLAT0)
      
      XE=ELON*DEG
      YE=ELAT*DEG
      
      ICL=1
      CON=-999.0      
      CALL LI2D(XE,YE,WA,HZM,'H',Z0,IC)       
      WA1(J)=Z0
      HMSK(J)=IC
      CALL LI2D(XE,YE,WA(1,2),HZM,'H',Z0,IC)
      WPP(J)=Z0      
      ENDDO
      
      WA1(NS)=WA(JMAX,1)
      WPP(NS)=WA(JMAX,2)      
      
C     END of COLUMN MASS at ALL SITES

      CALL INDEXS(JMAX,IMAX,IM,IP,JP)
      XP=XB+DEX*(IP-1)       
      YP=YB+DEY*(JP-1)
      ELON=XP*RADDEG
      ELAT=YP*RADDEG
      
      CALL FromETA(ELON,ELAT,GLON,GLAT,GLON0,GLAT0)       
      
      X(NS)=GLON*DEG
      Y(NS)=GLAT*DEG
      
C     END of MAXIMAL POINT COORDINATES      
C
C     Continue of Calculation
C
      DO J=1,IMJM
      C1=1000.*WA(J,2)/G
      DO K=1,LM
      DZ=(ZZ(J,K,1)-ZZ(J,K+1,1))
      WA(J,K)=C1*QQ(J,K,1)*HTM(J,K)*DETA(K)/DZ  ! mikro g/m^3
      ENDDO
      ENDDO 
      
C     END of DIMENSIONAL CONCENTRATION, micro g/m^3
      
      DO J=1,NS-1
               
      YJ=Y(J)
      GLAT=YJ*RADDEG           
      XI=X(J)
      GLON=XI*RADDEG 
            
      CALL ToETA(GLON,GLAT,
     *ELON,ELAT,GLON0,GLAT0)
      
      XE=ELON*DEG
      YE=ELAT*DEG
            
      DO K=1,LM
      CALL LI2D(XE,YE,WA(1,K),HZM,'H',Z0,IC)        !  mikro g/m^3  
      WCD(K,J)=Z0          
      ENDDO

      ENDDO
      
      DO K=1,KMAX
      WCD(K,NS)=WA(JMAX,K)          
      ENDDO
      
      DO J=1,NS
      CMAX=0.
      LMAX=1
            
      DO K=1,LM
      IF(WCD(K,J).GT.CMAX)THEN
      CMAX=WCD(K,J)
      LMAX=K
      ENDIF
      ENDDO
      
      WC(J)=CMAX
      WP(J)=0.01*(PT+AETA(LMAX)*WPP(J))    
      ENDDO      
            
      GAM=2.73
      
      WRITE(NU,100)GAM 

      WRITE(NU,*)
     *'             Site                    ',
     *'Lon   Lat    CM      OD      MC      LV'
    
      DO J=1,NS
      OD(J)=WA1(J)*GAM
      WRITE(NU,102)J,SITE(J),X(J),Y(J),WA1(J),OD(J),WC(J),WP(J)      ! ,HMSK(J)
      ENDDO
C23456789 123456789 123456789 123456789 123456789 123456789 123456789 12
      WRITE(NU,*)'    '      
      WRITE(NU,*)'        CM - Column mass, g/m^2'
      WRITE(NU,*)'        OD - Optical depth' 
      WRITE(NU,*)'        MX - Maximum concentration in column,',
     *                                  ' micro g/m^3' 
      WRITE(NU,*)'        LV - Level with maximum concentration, hPa'   
      WRITE(NU,*)'   '   
      WRITE(NU,*)'  Vertical Distribution (micro g/m^3)'

      DO K=1,LM
      WRITE(NU,103)K,(WCD(K,J),J=1,NS)      
      ENDDO                           

      
      RETURN    
100   FORMAT(10X,
     *      'Aerosol Mass (g/m^2) and Optical Depth for gamma=',
     *       F5.1,' m^2/g')

102   FORMAT(5X,I2,1X,A27,2(1X,F5.1),1X,2(1X,F7.4),2(1X,F7.1))
103   FORMAT(1X,I2,1X,10(1X,F6.1))              
      END
C& ATMOSM
      SUBROUTINE ATMOSM(NU,PS,ATM1,ATM2)
C
C     ATMOSPHERE MASS FOR TWO MOMENT (K,K-1)
C     IN CALCULATE POINT 
C
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)      
      PARAMETER (IMJM = IM*JM - JM/2)
      PARAMETER (KMAX = LM)
      PARAMETER (KKMAX = KMAX)      
      
      DIMENSION PS(IMJM,2)

      COMMON/ETAVEG/IHVEG(IMJM)
      COMMON/CBT1T2/HINTG(IMJM), HINN(IMJM,KMAX),
     *              HZM(IMJM,KMAX+1),BUFI(IMJM,2),      
     *              QFLUX(IMJM,2) ,EFIREA(IMJM,4),
     *              QSOUR(10),ISOUR(10),JSOUR(10),
     *              XSOUR(10),YSOUR(10),
     *              QFIL(IMJM,KKMAX,2),     
     *              IDATCP(3) 
      COMMON/ETAREA/CAREA(IMJM) 
      DATA G/9.81/
      
      SAM1=0.      !  total mass of air 
      SAM2=0.      !  total mass of air       
       
      DO J=1,IMJM
      C=HINTG(J)**CAREA(J)/G
      SAM1=SAM1+C*PS(J,1)          ! kg 
      SAM2=SAM2+C*PS(J,2)          ! kg       
      ENDDO

      ATM1=SAM1
      ATM2=SAM2
      
      WRITE(NU,*)'ATMOSPHERE MASS FOR TWO MOMENT',
     *             ' (K,K-1) IN CALCULATE POINT'
      WRITE(NU,*)'ATM(K-1)=',SAM1,' ATM(K)=',SAM2,' #/'
      RETURN
      END  
C& WTSCA
      SUBROUTINE WTSCA(Q,S,ACPR,DQDT,HTM,KMAX,DTT)
C
C           Wet scavenging (22/01/2001)
C
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)      
      PARAMETER (IMJM = IM*JM - JM/2)
      PARAMETER (LP1=LM+1, LM1=LM-1)
      PARAMETER (KKMAX=LM)
      DIMENSION Q(IMJM,KMAX),HTM(IMJM,KMAX),
     *          ACPR(IMJM,2),DQDT(IMJM,KMAX)
      COMMON/INNERP/IPN,INP(IMJM),IINP(IMJM),JINP(IMJM)
CX      
      COMMON/CBT1T2/HINTG(IMJM), HINN(IMJM,KKMAX),
     *              HZM(IMJM,KKMAX+1),BUFI(IMJM,2),      
     *              QFLUX(IMJM,2) ,EFIREA(IMJM,4),
     *              QSOUR(10),ISOUR(10),JSOUR(10),
     *              XSOUR(10),YSOUR(10),
     *              QFIL(IMJM,KKMAX,2),     
     *              IDATCP(3)            
      COMMON/ETAREA/CAREA(IMJM)
      DIMENSION R(IMJM,2)
      DATA SOLD/0.0046/   !     2)      
C     ----------------------------------------      
C      DATA S/0.00046/    !     1)
C      DATA S/0.0046/     !     2)
C     S=1 1/cm = 0.1 1/mm      
C     ----------------------------------------
C     Wet sedimentation rate R= S*(dP/dt)
C     (dP/dt) in mm/step;  1 step = 2880 sec
C     ACPR is accumulated precipitation in mm (kg/m**2)
C     S = 1/mm
C     ----------------------------------------
      S1=S
      DO5100 KK=1,IPN
      K=INP(KK)

      DQ=S1*(ACPR(K,2)-ACPR(K,1))      ! *1.E+6      !   into ppm
      DQOLD=SOLD*DTT*(ACPR(K,2)-ACPR(K,1))
      
      R(K,1)=DQ
      R(K,2)=DQOLD
      
      DO L=1,KMAX 
      
CX      DQ=S1*(ACPR(K,2)-ACPR(K,1))*HTM(K,L)*1.E+6
CX      DQOLD=SOLD*DTT*(ACPR(K,2)-ACPR(K,1))*HTM(K,L)
      
      QP=HTM(K,L)*Q(K,L)      
      QNEW=(QP-DQ)*HTM(K,L)
      QQ=AMAX1(QNEW,0.)
      DQ1=QP-QQ 
      DQDT(K,L)=HTM(K,L)*DQ1
      Q(K,L)=HTM(K,L)*QQ
      
      ENDDO
      
5100  CONTINUE

      CON =-999.0
      CALL AVSUR(R(1,1),HINTG,CON,RINT1,RAVE1)
      CALL AVSUR(R(1,2),HINTG,CON,RINT2,RAVE2)
      WRITE(66,*)' '
      WRITE(66,*)' Removal Rate (1/step)'
      WRITE(66,*)'S     =',S       ,' 1/mm'    
      WRITE(66,*)'Old S =',SOLD*DTT,'     '      
      WRITE(66,*)'RINT1=',RINT1,' RAVE1=',
     *                    RAVE1,' ppm/step (x 1.E+6)'
      WRITE(66,*)'RINT2=',RINT2,' RAVE2=',
     *                    RAVE2,' ppm/step (old)'
      
      CALL EXTR2M(66,
     * 'SCAVENING PER STEP (ppm/step)',
     * DQDT,HTM,HINTG,IMJM,KMAX,1)
      
C      CALL EXTR2M(66,'ETADOT', WW(1,1,2),HTM(1,KTOP),HINTG,
C     *          IMJM,KMAX,1)                 
      RETURN
      END
C
C               One fire subroutines
C          
C& TRANXY
      SUBROUTINE TRANXY(HX,X1,X2,X3,HY,Y1,Y2,Y3,FLON,FLAT)
      FLON=HX*(X1+X2/60.+X3/3600.)
      FLAT=HY*(Y1+Y2/60.+Y3/3600.)
      RETURN
      END
C& ONEFIR
      SUBROUTINE ONEFIR(EFIREA,ARFKM2,FLON,FLAT)
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)      
      PARAMETER (IMJM = IM*JM - JM/2)
      PARAMETER (LENT=1000)
      
      DIMENSION EFIREA(IMJM,4)
      

      RADDEG=3.14159265/180.
      DEG=180./3.14159265             
      
      CALL CGLOB(GLON0,GLAT0)

      WRITE(66,*)'ONEFIR: FLON, FLAT =',FLON,FLAT
      WRITE(66,*)'ONEFIR: GLON0,GLAT0=',GLON0,GLAT0
      
            
      FLON1=FLON*RADDEG 
      FLAT1=FLAT*RADDEG
      CALL ToETA(FLON1,FLAT1,ELON1,ELAT1,GLON0,GLAT0)
      ELON=ELON1*DEG 
      ELAT=ELAT1*DEG 
      IPOINT=INDEXC(ELON,ELAT)  
      WRITE(66,*)'ONEFIR:      IPOINT=',IPOINT
          

        IF(IPOINT.LT.1.OR.IPOINT.GT.IMJM)THEN
        WRITE(6,*)'ERROR:'
        WRITE(6,*)'FLON=',FLON,' FLAT=',FLAT
        WRITE(6,*)'ELON=',ELON,' ELAT=',ELAT
        WRITE(6,*)'IPOINT=',IPOINT
        STOP 'POINT ERROR'        
        ENDIF
        

      EFIREA(IPOINT,1)=EFIREA(IPOINT,1)+ARFKM2       
      
      RETURN      
      END
C& CGLOB
      SUBROUTINE CGLOB(GLON0,GLAT0)           
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)
      PARAMETER (LP1=LM+1)            
      PARAMETER (IMJM = IM*JM - JM/2)
      PARAMETER (JM2=JM/2+1,IMT=IM*2-1)      
      PARAMETER (NPOINT=25)
      PARAMETER (IMAX=51, JMAX=51)      
C
      COMMON/ETAGRD/
     *      XB,DEX,YB,DEY,
     *      COSETA(JM),XETA(IMT),YETA(JM),
     *      ITMP(NPOINT), JTMP(NPOINT), KTMP(NPOINT)
C----------------------------------------------------
                              C O M M O N /MAPOT/
     & LSL,IXM,IYM
     &,ISHDE(49)
     &,R,PT,TSPH
     &,WBD,SBD,TLM0D,TPH0D,DLMD,DPHD,CMLD,DP30,X1P,Y1P
     &,DISLP,Z0SLP,SPL(LSM),ALSL(LSM)
     &,DETA(LM),AETA(LM),ETA(LP1),TSHDE(49)
     &,ERLAM0,CPHI0,SPHI0
     &,KUNIT,LUNIT,MUNIT,IUNIT,NUNIT
C-----------------------------------------------------      
      DATA ERAD /6370000./
C      
       RADDEG=3.14159265/180.
       DEG=180./3.14159265
C
       DPHI=DPHD*RADDEG   !   in radian 
       DLAM=DPHD*RADDEG   !   in radian 
       
       GLON0=TLM0D*RADDEG
       GLAT0=TPH0D*RADDEG
       RETURN
       END 
C& FLUXON
      SUBROUTINE FLUXON(ICLF1,QFLUX,QFLX,EFIREA)
C
C     SURFACE FLUX
C
      PARAMETER (IM=111,JM=221,LM=38,LSM=19)      
      PARAMETER (IMJM = IM*JM - JM/2)
      PARAMETER (LP1=LM+1)
      PARAMETER (JAM=6+2*(JM-10))      
      PARAMETER (IMAX=2*IM-1  , JMAX=JM)       
      PARAMETER (KKMAX=LM)
      PARAMETER (KMAX=KKMAX)
      PARAMETER (KMIN = LM-KKMAX)
      
                 
      DIMENSION QFLUX(IMJM,2),EFIREA(IMJM,4)
C-----------------------------------------------------------------------
                             C O M M O N /LOOPS/
     & KHLA  (JAM),KHHA  (JAM),KVLA  (JAM),KVHA  (JAM)
C
     &,KHL2  (JM),KHH2  (JM),KVL2  (JM),KVH2  (JM)
C
     &,LMH   (IMJM),LMV   (IMJM)     
C------------------------------------------------------------------------
                              C O M M O N /MAPOT/
     & LSL,IXM,IYM
     &,ISHDE(49)
     &,R,PT,TSPH
     &,WBD,SBD,TLM0D,TPH0D,DLMD,DPHD,CMLD,DP30,X1P,Y1P
     &,DISLP,Z0SLP,SPL(LSM),ALSL(LSM)
     &,DETA(LM),AETA(LM),ETA(LP1),TSHDE(49)
     &,ERLAM0,CPHI0,SPHI0
     &,KUNIT,LUNIT,MUNIT,IUNIT,NUNIT
C-----------------------------------------------------------------------
                             C O M M O N /VRBLS/
     & PD(IMJM),RES(IMJM),FIS(IMJM)
C
     &,U(IMJM,LM),V(IMJM,LM),T(IMJM,LM),Q(IMJM,LM)
C-----------------------------------------------------------------------
      COMMON/ETATRN/
     * QQL(IMJM,KKMAX),QQH(IMJM,KKMAX),      
     * PS(IMJM,2), QQ(IMJM,KKMAX,2),  
     * UU(IMJM,KKMAX,2), VV(IMJM,KKMAX,2),
     * WW(IMJM,KKMAX,2), RO(IMJM,KKMAX,2),     
     * TT(IMJM,KKMAX,2), ZZ(IMJM,KKMAX+1,2),
     * ALFA(IMJM,KKMAX,2),
     * DIFQ(IMJM,KKMAX),
     * AKH(IMJM,KKMAX,2),AKHS(IMJM,2),
     * TOPG(IMJM)
     *,PREC(IMJM,2),ACPR(IMJM,2)
C-----------------------------------------------------------------------
C      
      COMMON/ETAREA/CAREA(IMJM)
C
      DATA RGAS/287.04/    
C-----------------------------------------------------------------------          
      IF(ICLF1.GT.1)THEN
      WRITE(66,*)'Error in FLUXON: ICLF1=',ICLF1
      STOP
      ENDIF
C
C                 QFLUX(*,1)
C      
C
C     EFIREA(J) in km**2
C     CAREA(J)  in m**2
C
      C1=QFLX*1.E6  ! For converting EFIREA(J) to m**2
      C2=C1*1.E6    ! For converting to ppm          
      DO J=1,IMJM
      QFLUX(J,1)=ICLF1*QFLUX(J,1)+C2*EFIREA(J,1)/CAREA(J)
      EFIREA(J,2)=0.
      EFIREA(J,3)=0.
      EFIREA(J,4)=0.
      ENDDO      
C
C                 QFLUX(*,2)
C
       DO10 J=1,IMJM
       
       C=PD(J)*RES(J)
       K=LMH(J)
       IF(K.EQ.0)THEN
       QFLUX(J,2)=0.
       QFLUX(J,1)=0.
       ELSE
       PSUR=PT+AETA(K+KMIN)*C
       TEM=TT(J,K,2)
       ROP=PSUR/(RGAS*TEM)
       DZ=ZZ(J,K,2)-ZZ(J,K+1,2)
       QFLUX(J,2)=QFLUX(J,1)/(ROP*DZ)
       ENDIF
       
 10    CONTINUE
C      
      RETURN
      END
C& RDBUFx
      SUBROUTINE RDBUFx(NU,ISTEP,IREC,BUFI,FIRE,IMJM)
      DIMENSION BUFI(IMJM,2),FIRE(IMJM)
      DIMENSION INX(30), ALPHA(4)
      DATA INX/1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,
     *         3,3,3,3,3,3,3,4,4,4,4,4,4,4,4/
      DATA ALPHA/0.47,0.53,0.56,0.44/
         
      IDAY   = (ISTEP-1)/30+1      !             day
      TIME   = (IDAY-1)*24.*60.    !      min    time
      
      ISTED  = MOD(ISTEP-1,30)+1   !             inner day step
      TIMED  = ISTED*48.           !      min    inner time
      WRITE(66,*)'RDBUF: ISTEP=',ISTEP,' step'     
      WRITE(66,*)'RDBUF: IDAY =',IDAY, ' day'
      WRITE(66,*)'RDBUF: TIME =',TIME, ' min'
      WRITE(66,*)'RDBUF: ISTED=',ISTED,' inner day, step' 
      WRITE(66,*)'RDBUF: TIMED=',TIMED,' inner time, min'
             
          IF(ISTED.EQ.1.OR.ISTED.EQ.16)THEN       !  new
!                                   15            !  old   	 
	 READ(NU,rec=IREC)(BUFI(J,1),J=1,IMJM)       
      WRITE(66,*)'RDBUFx: read IREC=',IREC,' to BUFI(J,1)'	 	            
      WRITE(66,*)
     *'FIRE AREA OF IDAY=',IDAY,' HOURE=',TIME,      
     * ' IREC=',IREC,' /#' 
	 IREC=IREC+1
	 
	 IF(IREC.EQ.3)IREC=1    !  new   Next day is percisten (second)
	                      	 	 
	 ENDIF
	 
	 
             IA=INX(ISTED)        !  new
             ALP=ALPHA(IA)        !  new
             	 
             DO J=1,IMJM
             FIRE(J)=ALP*BUFI(J,1)  !  new            
             ENDDO		  	 	                        	  
      RETURN   
      END
C& LOADMS       
      SUBROUTINE LOADMS(NU,QFLUX,EFIREA,IMJM,DT)
      DIMENSION  QFLUX(IMJM),EFIREA(IMJM)
      REAL SLOAD,S
      S=0.
      AREA=0.      
      DO J=1,IMJM
      S=S+QFLUX(J)*EFIREA(J)
      AREA=AREA+EFIREA(J)
      ENDDO
      
      
      SLOAD=S*1.E6*DT  !    (1.E6 transforms km**2 to m**2)
      
      WRITE(NU,100)SLOAD,DT,AREA
      RETURN
C23456
100   FORMAT(' Loaded mass=',E12.5,' kg, during ',F7.1,' sec, from ',
     & F10.5,' km**2')
      END
C& EXDIAG                                                               
      SUBROUTINE EXDIAG(NU,NAME,A,IMJM,KMAX) 
      PARAMETER (LM=38) 
      DIMENSION A(IMJM,KMAX) 
      CHARACTER*(*) NAME 
      DIMENSION B1(LM),B2(5),K1(LM),K2(5) 
      CHARACTER*6 LEV(5) 
      DATA LEV/'    1',' 2-13','14-25','26-37','   38'/ 
C                   S T A R T                                           
      WRITE(66,*)'EXDIAG for 1d ',NAME 
      IF(KMAX.GT.LM)THEN
      WRITE(66,*)'Error in EXDIAG: KMAX=',KMAX,' > 38'
      STOP
      ENDIF
                                                                        
      DO K=1,KMAX                 !  DO 1000 
      B1(K)=0. 
      K1d=0 
      S1d=0. 
                                                                        
      DO J=1,IMJM                 !  DO 1001
      S1d=S1d+A(J,K) 
      K1d=K1d+1 
      END DO                      !   1001 END DO
                                                                        
      S1d=S1d/FLOAT(K1d) 
      B1(K)=S1d 
      K1(K)=K1d 
      END DO                      ! 1000 END DO  
                                                                        
      IF(KMAX.EQ.1)THEN 
      WRITE(66,100)1,' 1d= ',K1(1),B1(1),' IMJM=',IMJM 
      RETURN 
      ENDIF 
                                                                        
      IF(KMAX.EQ.38)THEN 
         DO K=2,4 
         B2(K)=0. 
         K2(K)=0 
         ENDDO 
                                                                        
      C1=1./12. 
      B2(1)=B1(1) 
      B2(5)=B1(38)
      K2(1)=K1(1) 
      K2(5)=K1(38)        
                                                                        
      DO K=1,12 
      B2(2)=B2(2)+B1(1+K)*C1 
      B2(3)=B2(3)+B1(13+K)*C1 
      B2(4)=B2(4)+B1(25+K)*C1 
      K2(2)=K2(2)+K1(1+K) 
      K2(3)=K2(3)+K1(13+K) 
      K2(4)=K2(4)+K1(25+K) 
      ENDDO 
                                                                        
         DO K=2,4 
         K2(K)=K2(K)/12 
         ENDDO 
                                                                        
      DO K=1,5 
      WRITE(66,101)LEV(K),' 1d= ',K2(K),B2(K),' IMJM=',IMJM 
      ENDDO 
                                                                        
      RETURN 
      ENDIF 
                                                                        
      IF(KMAX.NE.1.AND.KMAX.NE.38)THEN 
      DO K=1,KMAX 
      WRITE(66,100)K,' 1d= ',K1(K),B1(K),' IMJM=',IMJM 
      ENDDO 
      ENDIF 
                                                                        
      RETURN 
100   FORMAT(I3,4X,A5,I6,1X,E11.4,A6,3(1X,I5)) 
101   FORMAT(A6,1X,A5,I6,1X,E11.4,A6,3(1X,I5)) 
      END  
                       
