      PROGRAM DYNNH3EMIS
C
C **************************************************************************
C * THIS PROGRAM IS THE SECOND VERSION OF THE DYNAMICAL AMMONIA            *
C * AMMONIA PARAMETERISATION                                               *
C *                                                                        *
C * INPUT: GRIDDED HOURLY METEOROLOGY: TEMP WIND SPEED
C *                                                                        *
C * DOCUMENTATION: "IMPLEMENTING A DYNAMICAL AMMONIA PARAMETERIZATION      *
C *                 IN THE LARGE-SCALE AIR POLLUTION MODEL ACDEP"          *
C *                 C. A. SKJOETH, O. HERTEL, S. GYLDENKAERNE AND          *
C *                 T. ELLERMANN, JGR-ATMOSPHERE, VOL 109, D06306          *
C *                 doi:10.1029/2003JD003895,2004                          *
C *                                                                        *
C * DOCUMENTATION: "A DYNAMICAL AMMONIA EMISSION PARAMETERIZATION          *
C *                 FOR USE IN AIR POLLUTION MODELS,                       *
C *                 S. GYLDENKAERNE, C.A. SKJOETH, O. HERTEL AND           *
C *                 T.ELLERMANN, JGR-ATMOSPHERE, VOL 110, D07108           *
C *                 doi:10.1029/2004JD005459,2005                          *
C *                                                                        *
C *                "FOOTPRINTS ON AMMONIA FROM ENVIRONMENTAL REGULATIONS,  *
C *                 C. A. SKJOETH, T.ELLERMANN, O.HERTEL, S.GYLDENKAERNE   *
C *                 AND M. H. MIKKELSEN, J. AIR & WASTE MANAGEMENT ASSOC.  *
C *                 58:1158-1165, 2008.                                    *
C *                 doi:10.3155/1047-3289.58.9.1158                        *
C *                                                                        *
C *                "SPATIAL AND TEMPORAL VARIATIONS IN AMMONIA EMSSIONS -  *
C *                 A FREELY ACCESSIBLE MODEL CODE FOR EUROPE              *
C *                 C.A.SKJOETH ET AL (2010), SUBM TO                      *
c *                 ATMOPHERIC CHEMISTRY AND PHYSCIS, NOV 20, 2010         *
C *                                                                        *
C * COMPILATION: 'pgf77 -byteswapio -o nh3emis_v2.x nh3emis.ver2.f         *
C *                                                                        *
C * PREPARED FOR THE DEHM AND EMEP MODELS MARCH 2010                       *
C * C.A.SKJOETH MARCH 2010                                                 *
C **************************************************************************
C
      IMPLICIT NONE
C
      INTEGER I, J, IMONTH(12), NDAYS, JJ,
     &        ISTARTYEAR, LL, II, III, IGRIDNR
      INTEGER ID, IDATE, IDAYTAELLER, 
     &        IIDAYS, KK,L,M,
     &        KKK, NHOURS, NXNY
      PARAMETER (IIDAYS=366*24,NXNY=96*96)
      INTEGER DAYGRIDTEMP(8,NXNY)
      REAL RTEMP, T2MEAN(366), K(IIDAYS), 
     &     VH10GRID(NXNY), TEMPGRID(NXNY), TEMP_DAY_MEAN(NXNY),
     &     T2MEAN_GRID(366,NXNY),  
     &     SUMT2GRID(NXNY), U10GRID(NXNY), V10GRID(NXNY)
      REAL RTEMPa(NXNY), SUMISOLATED(NXNY), SUMNONISOLATED(NXNY),
     &     SUMSTORAGE(NXNY), W10GRID(NXNY), W10HELPGRID(NXNY),
     &     SFKTGRID8(NXNY), SFKTGRID9(NXNY), SFKTGRID10(NXNY),
     &     SFKTGRID11(NXNY), SFKTGRID11a(NXNY), SFKTGRID12(NXNY),
     &     SFKTGRID13(NXNY), SFKTGRID14(NXNY), SFKTGRID4(NXNY),
     &     SFKTGRID5(NXNY), SFKTGRID6(NXNY), SFKTGRID7(NXNY)
C
      REAL FKT1GRID(NXNY), FKT2GRID(NXNY), FKT3GRID(NXNY),
     &     FKT8GRID(NXNY),  FKT9GRID(NXNY), FKT10GRID(NXNY),
     &     FKT11GRID(NXNY), FKT11aGRID(NXNY), FKT14GRID(NXNY),
     &     FKT4GRID(NXNY), FKT5GRID(NXNY), FKT6GRID(NXNY),
     &     FKT7GRID(NXNY), FKT12GRID(NXNY), FKT13GRID(NXNY),
     &     FKT15GRID(NXNY), FKT16GRID(NXNY),FKTALLGRID(NXNY)
C
      REAL GRAZING(NXNY), GRASSGROWTH(NXNY), SUMGRAZING(NXNY), 
     &     IGRASSGROWTH(7*24,NXNY), GRAZING_GRID(NXNY)

      CHARACTER*28 FILENAME
      CHARACTER*4 CYEAR
      CHARACTER*2 CMONTH
      CHARACTER*6 TYPE
      LOGICAL VARIATION, DUMB
      PARAMETER (VARIATION=.FALSE.) !TRUE EMIS VARIATION, FALSE ACTUAL EMISSIONS
      PARAMETER (DUMB=.FALSE.) !TRUE DUMB HOURLY GRIDDED FILES, FALSE NO DUMB - SAVING TIME

      INTEGER JJJ
      REAL RHELP, RHELP2(17)
      CHARACTER*4 EMISYEAR*4
      REAL EAST, NORTH, LAT, LON, EMIS(16), 
     &     GEMEMIS(16,NXNY), TIME1
      REAL C1, D1, PI,SQRTTWOPI
C
      REAL TSTALD, VENTFACT
C
      DATA IMONTH/ 31,28,31,30,31,30,31,31,30,31,30,31/
      PI=4.*ATAN(1.)
      SQRTTWOPI=(PI*2.)**0.5
      READ(*,*) ISTARTYEAR
      WRITE(CYEAR,'I4') ISTARTYEAR
      IF (ISTARTYEAR.EQ.1980.OR.ISTARTYEAR.EQ.1984.OR.
     &    ISTARTYEAR.EQ.1988.OR.ISTARTYEAR.EQ.1992.OR.
     &    ISTARTYEAR.EQ.1996.OR.ISTARTYEAR.EQ.2000.OR.
     &    ISTARTYEAR.EQ.2004.OR.ISTARTYEAR.EQ.2008.OR.
     &    ISTARTYEAR.EQ.2012.OR.ISTARTYEAR.EQ.2016) THEN
          IMONTH(2)=29
          NDAYS=366
      ELSE
          NDAYS=365
      ENDIF
      NHOURS=NDAYS*24
      WRITE(EMISYEAR,'(I4)') ISTARTYEAR
C     
C *******************************************************************
C * CONTROLLING THE EMISSION MODEL:                                 *
C * EITHER NORMALIZED VARIATIONS                                    *
C * OR ACTUAL EMISSSIONS BASED ON AN INPUT FILE WITH                *
C * GRIDDED VALUES FROM EACH SECTION(FUNKTIONS)                     *
C *******************************************************************
C
      IF (VARIATION) THEN
      WRITE(*,*) 'EMISVAR'
      DO II=1,NXNY
      DO I=1,16
         GEMEMIS(I,II)=1.
      ENDDO
      ENDDO
      ELSE
      WRITE(*,*) 'ACTUAL EMISSIONS'

      OPEN(31,FILE='Sector_Emis.txt',TYPE='UNKNOWN')
      READ(31,*)
      DO II=1,NXNY
      READ(31,*) IGRIDNR,(RHELP2(I),I=1,17)
      GEMEMIS(01,IGRIDNR)=RHELP2(01)
      GEMEMIS(02,IGRIDNR)=RHELP2(02)
      GEMEMIS(03,IGRIDNR)=RHELP2(03)
      GEMEMIS(04,IGRIDNR)=RHELP2(04)
      GEMEMIS(05,IGRIDNR)=RHELP2(05)
      GEMEMIS(06,IGRIDNR)=RHELP2(06)
      GEMEMIS(07,IGRIDNR)=RHELP2(07)
      GEMEMIS(08,IGRIDNR)=RHELP2(08)
      GEMEMIS(09,IGRIDNR)=RHELP2(09)
      GEMEMIS(10,IGRIDNR)=RHELP2(10)
      GEMEMIS(11,IGRIDNR)=RHELP2(11)
      GEMEMIS(12,IGRIDNR)=RHELP2(12)
      GEMEMIS(13,IGRIDNR)=RHELP2(13)+RHELP2(14)
      GEMEMIS(14,IGRIDNR)=RHELP2(15)
      GEMEMIS(15,IGRIDNR)=RHELP2(16)
      GEMEMIS(16,IGRIDNR)=RHELP2(17)
      ENDDO
      ENDIF
C
      I=0
      JJ=0
      OPEN(13,FILE='emisdata.test.2009.dat')
      OPEN(14,FILE='out.test.dat')
      OPEN(15,FILE='Grazing_days_gridded16km.prn')
      DO II=1,NXNY
         READ(15,*) IGRIDNR,RHELP
         IF (RHELP.LE.0) THEN
             GRAZING_GRID(IGRIDNR)=183.
         ELSE
             GRAZING_GRID(IGRIDNR)=RHELP
         ENDIF
      ENDDO
      CLOSE(15)
C METEOROLOGICAL INPUT FILES
      OPEN(22,FILE=
     &'/nfs2data03/cas/dyn.pollen.emis/input/T2FIELD.'//CYEAR//'.N2',
     &TYPE='OLD',FORM='UNFORMATTED')
C
  307 CONTINUE
      DO L=1,365
      DO II=1,NXNY
         TEMP_DAY_MEAN(II)=0.
      ENDDO
      DO KKK=1,24
      JJ=JJ+1
      READ(22) TEMPGRID
      DO II=1,NXNY
         TEMPGRID(II)=TEMPGRID(II)-273.14
         TEMP_DAY_MEAN(II)= TEMPGRID(II)/24+TEMP_DAY_MEAN(II)
      ENDDO
C
C *******************
C 
      ENDDO !KKK
      I=I+1 !DAYNUMBER
      DO II=1,NXNY
         T2MEAN_GRID(I,II)=TEMP_DAY_MEAN(II)
      ENDDO
      ENDDO !L
  308 CONTINUE
      CLOSE(12)
C
      IDAYTAELLER=IMONTH(1)+IMONTH(2)
      DO II=1,NXNY
      SUMT2GRID(II)=0.
      DAYGRIDTEMP(1,II)=0.
      DAYGRIDTEMP(2,II)=0.
      DAYGRIDTEMP(3,II)=0.
      DAYGRIDTEMP(4,II)=0.
      DAYGRIDTEMP(5,II)=0.
      DAYGRIDTEMP(6,II)=0.
      DAYGRIDTEMP(7,II)=0.
      DAYGRIDTEMP(8,II)=0.
      ENDDO
      DO I=IMONTH(1)+IMONTH(2),NDAYS
         IDAYTAELLER=IDAYTAELLER+1
         DO II=1,NXNY
            SUMT2GRID(II)=SUMT2GRID(II)+T2MEAN_GRID(IDAYTAELLER,II)
            IF (SUMT2GRID(II).GT.090.AND.DAYGRIDTEMP(1,II).EQ.0) THEN
             DAYGRIDTEMP(1,II)=IDAYTAELLER
            ENDIF
            IF (SUMT2GRID(II).GT.200.AND.DAYGRIDTEMP(2,II).EQ.0) THEN
             DAYGRIDTEMP(2,II)=IDAYTAELLER
            ENDIF
            IF (SUMT2GRID(II).GT.260.AND.DAYGRIDTEMP(3,II).EQ.0) THEN
             DAYGRIDTEMP(3,II)=IDAYTAELLER
            ENDIF
            IF (SUMT2GRID(II).GT.1700.AND.DAYGRIDTEMP(4,II).EQ.0) THEN
             DAYGRIDTEMP(4,II)=IDAYTAELLER
            ENDIF
            IF (SUMT2GRID(II).GT.710.AND.DAYGRIDTEMP(5,II).EQ.0) THEN
             DAYGRIDTEMP(5,II)=IDAYTAELLER
            ENDIF
            IF (SUMT2GRID(II).GT.950.AND.DAYGRIDTEMP(6,II).EQ.0) THEN
             DAYGRIDTEMP(6,II)=IDAYTAELLER
            ENDIF
            IF (SUMT2GRID(II).GT.1400.AND.DAYGRIDTEMP(7,II).EQ.0) THEN
             DAYGRIDTEMP(7,II)=IDAYTAELLER
            ENDIF
            IF (SUMT2GRID(II).GT.2040.AND.DAYGRIDTEMP(8,II).EQ.0) THEN
             DAYGRIDTEMP(8,II)=IDAYTAELLER
            ENDIF


         ENDDO
      ENDDO
      WRITE(*,*) 'WRITING GAUSS PEAKS'
      OPEN(11,FILE='test.dat')
      DO II=1,NXNY
      WRITE(11,'(9I6)') II,DAYGRIDTEMP(1,II),DAYGRIDTEMP(2,II),
     & DAYGRIDTEMP(3,II),DAYGRIDTEMP(4,II),DAYGRIDTEMP(5,II),
     & DAYGRIDTEMP(6,II),DAYGRIDTEMP(7,II),DAYGRIDTEMP(8,II)
      ENDDO
      CLOSE(11)
      CLOSE(21)
      CLOSE(22)
      OPEN(21,FILE=
     &'/nfs2data03/cas/dyn.pollen.emis/input/W10FIELD.'//CYEAR//'.N2',
     &TYPE='OLD',FORM='UNFORMATTED')
      OPEN(22,FILE=
     &'/nfs2data03/cas/dyn.pollen.emis/input/T2FIELD.'//CYEAR//'.N2',
     &TYPE='OLD',FORM='UNFORMATTED')
      WRITE(*,*)'PREPROCESSING EMISSION POTENTIALS'
C
C *******************************
C * INITIALIZATION OF VARIABLES *
C *******************************
C 
      DO II=1,NXNY
         GRAZING(II)=0
         GRASSGROWTH(II)=0
         SUMISOLATED(II)=0.
         SUMNONISOLATED(II)=0.
         SUMSTORAGE(II)=0.
         SFKTGRID4(II)=0.
         SFKTGRID5(II)=0.
         SFKTGRID6(II)=0.
         SFKTGRID7(II)=0.
         SFKTGRID8(II)=0.
         SFKTGRID9(II)=0.
         SFKTGRID10(II)=0.
         SFKTGRID11(II)=0.
         SFKTGRID11a(II)=0.
         SFKTGRID12(II)=0.
         SFKTGRID13(II)=0.
         SFKTGRID14(II)=0.
      DO JJJ=1,7*24
      IGRASSGROWTH(JJJ,II)=0.
      ENDDO

      ENDDO
C
      TIME1=0
      DO L=1,365
      DO KKK=1,24
      TIME1=TIME1+1
      READ(21) U10GRID
      READ(21) V10GRID
      READ(22) TEMPGRID
      DO II=1,NXNY
         VH10GRID(II)=SQRT(U10GRID(II)**2.+V10GRID(II)**2.)
         TEMPGRID(II)=TEMPGRID(II)-273.14
C
C *******************************************
C * EMISSION POTENTIAL FOR ISOLATED STABLES *
C *******************************************
C
      IF (TEMPGRID(II).LT.0) THEN
      VENTFACT=0.2
      TSTALD=MAX(0.,(18.0+0.5*(TEMPGRID(II)-0.)))
      ELSE IF (TEMPGRID(II).GE.0.AND.TEMPGRID(II).LT.12.5) THEN
      TSTALD=18.
      VENTFACT=(0.2+(0.38-0.2)*(MAX(0.,TEMPGRID(II))-0.)/(12.5-0.))
      ELSE
      VENTFACT=0.38
      TSTALD=18.+(TEMPGRID(II)-12.5)*0.77
      ENDIF
      SUMISOLATED(II)=SUMISOLATED(II)+
     & (VENTFACT**0.26)*TSTALD**0.89
C
C ***************************************
C * EMISSION POTENTIAL FOR OPEN STABLES *
C ***************************************
C
      IF (TEMPGRID(II).LT.1.) THEN
      TSTALD=4.
      ELSE
      TSTALD=TEMPGRID(II)+3.
      ENDIF
      VENTFACT=0.214+
     & 0.014*SIN((243.*TIME1+(REAL(TIME1)))*2.*PI/(24.*REAL(NDAYS)))
      SUMNONISOLATED(II)=SUMNONISOLATED(II)+
     &                   (VENTFACT**0.26)*(TSTALD**0.89)
C
C **********************************************
C * EMISSION POTENTIAL FOR STORAGE FASCILITIES *
C **********************************************
C
      IF (TEMPGRID(II).LT.1.) THEN
      TSTALD=1.
      ELSE
      TSTALD=TEMPGRID(II)
      ENDIF
      W10GRID(II)=VH10GRID(II)
      SUMSTORAGE(II)=SUMSTORAGE(II)+(W10GRID(II)**0.26)*TSTALD**0.89
C
      W10HELPGRID(II)=EXP(0.0419*W10GRID(II))
C
C ******************************************************
C * EMISSION POTENTIAL FOR MANURE TYPE 1: EARLY SPRING *
C ******************************************************
C
      C1=(DAYGRIDTEMP(1,II))*24.
      D1=5.*24.
      D1=9.*24.
      SFKTGRID8(II)=EXP(0.0223*TEMPGRID(II))*(1./(D1*SQRTTWOPI))*
     &      W10HELPGRID(II)*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)+SFKTGRID8(II)
C
C *****************************************************
C * EMISSION POTENTIAL FOR MANURE TYPE 2: LATE SPRING *
C *****************************************************
C
      C1=(DAYGRIDTEMP(3,II))*24.
      D1=5.*24.
      D1=9.*24.
      SFKTGRID9(II)=EXP(0.0223*TEMPGRID(II))*(1./(D1*SQRTTWOPI))*
     &      W10HELPGRID(II)*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)+SFKTGRID9(II)
C 
C ************************************************
C * EMISSION POTENTIAL FOR MANURE TYPE 3: SUMMER *
C ************************************************
C      
      C1=(DAYGRIDTEMP(7,II)) *24.
      D1=16.*24.
      SFKTGRID10(II)=EXP(0.0223*TEMPGRID(II))*(1./(D1*SQRTTWOPI))*
     &      W10HELPGRID(II)*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)+SFKTGRID10(II)
C
C ************************************************
C * EMISSION POTENTIAL FOR MANURE TYPE 4: AUTUMN *
C ************************************************
C
      C1=(270.)*24.
      D1=16.*24.
      SFKTGRID11(II)=EXP(0.0223*TEMPGRID(II))*(1./(D1*SQRTTWOPI))*
     &      W10HELPGRID(II)*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)+SFKTGRID11(II)
C
C ******************************************************
C * EMISSION POTENTIAL FOR MANURE TYPE 4a: EMPTY TANKS *
C ******************************************************
C
      C1=(270.)*24.
      D1=9.*24.
      SFKTGRID11a(II)=EXP(0.0223*TEMPGRID(II))*(1./(D1*SQRTTWOPI))*
     &      W10HELPGRID(II)*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)+SFKTGRID11a(II)
C
C *****************************************************
C * EMISSION POTENTIAL FOR MINERAL FERTILIZER: SPRING *
C *****************************************************
C
      C1=(DAYGRIDTEMP(2,II))*24.
      D1=9.*24.
      SFKTGRID12(II)=EXP(0.0223*TEMPGRID(II))*(1./(D1*SQRTTWOPI))*
     &      W10HELPGRID(II)*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)+SFKTGRID12(II)
C
C *****************************************************
C * EMISSION POTENTIAL FOR MINERAL FERTILIZER: SUMMER *
C *****************************************************
C
      C1=(DAYGRIDTEMP(7,II))*24.
      D1=16.*24.
      SFKTGRID13(II)=EXP(0.0223*TEMPGRID(II))*(1./(D1*SQRTTWOPI))*
     &      W10HELPGRID(II)*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)+SFKTGRID13(II)
C ******************************************************************************
C * EMISSION POTENTIAL FOR GRAZING CATTLE: ASSUMED TO FOLLOW THE GRASS PATTERN *
C ******************************************************************************
C
      C1=DAYGRIDTEMP(7,II)*24.
      D1=60.*24.
      D1=(GRAZING_GRID(II)/2)*24 !ASSUMING 95% COVERS THE TABULAR VALUE
      SFKTGRID14(II)=(1./(D1*SQRTTWOPI))*
     &      EXP(0.0223*TEMPGRID(II))*W10HELPGRID(II)*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)+SFKTGRID14(II)

cc      IF (TEMPGRID(II).LT.1.) THEN
cc      TSTALD=1.
cc     ELSE
cc      TSTALD=TEMPGRID(II)
cc      ENDIF
cc      IF (TEMPGRID(II).LT.0.) THEN
cc      GRAZING(II)=0
cc      GRASSGROWTH(II)=0
cc      ELSE
cc      GRASSGROWTH(II)=GRASSGROWTH(II)+TEMPGRID(II)-IGRASSGROWTH(1,II)
cc      ENDIF
C
C ******************************************************************
C * ACCUMULATES GRASSGROWTH LAST 7 DAYS - THE GRAZING REQUIREMENT *
C ******************************************************************
C
cc      DO JJJ=2,7*24
cc      IGRASSGROWTH(JJJ-1,II)=IGRASSGROWTH(JJJ,II)
cc      ENDDO
cc      IGRASSGROWTH(7*24,II)=TEMPGRID(II)
cc
cc      IF (GRASSGROWTH(II).GE.500.AND.
cc     &   (TIME1.GT.15*24.AND.TIME1.LT.350*24)) THEN
cc
cc      GRAZING(II)=1
cc      SUMGRAZING(II)=SUMGRAZING(II)+1
cc      ELSE
cc      GRAZING(II)=0
cc     ENDIF
cc      SFKTGRID14(II)=GRAZING(II)*(W10GRID(II)**0.26)*TSTALD**0.89+
cc     &               +SFKTGRID14(II)
C
C ***************************************
C * EMISSION POTENTIAL FOR WINTER CROPS *
C ***************************************
C
      C1=DAYGRIDTEMP(5,II)*24.
      D1=39.*24.
      SFKTGRID4(II)=(1./(D1*SQRTTWOPI))*
c      SFKTGRID4(II)=(1./(D1*SQRTTWOPI))*GEMEMIS(1,II)*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)+SFKTGRID4(II)
C
C ***************************************
C * EMISSION POTENTIAL FOR SPRING CROPS *
C ***************************************
C
      C1=DAYGRIDTEMP(6,II)*24.
      D1=25.*24.
      SFKTGRID5(II)=(1./(D1*SQRTTWOPI))*
c      SFKTGRID5(II)=(1./(D1*SQRTTWOPI))*GEMEMIS(2,II)*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)+SFKTGRID5(II)
C 
C *********************************************
C * EMISSION POTENTIAL FOR SPRING SUGAR BEETS *
C *********************************************
C
      C1=DAYGRIDTEMP(8,II)*24.
      D1=45.*24.
      SFKTGRID6(II)=(1./(D1*SQRTTWOPI))*
c      SFKTGRID6(II)=(1./(D1*SQRTTWOPI))*GEMEMIS(3,II)*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)+SFKTGRID6(II)
C 
C ********************************************
C * EMISSION POTENTIAL FOR GRASS IN ROTATION *
C ********************************************
C
      C1=DAYGRIDTEMP(7,II)*24.
      D1=60.*24.
      SFKTGRID7(II)=(1./(D1*SQRTTWOPI))*
c      SFKTGRID7(II)=(1./(D1*SQRTTWOPI))*GEMEMIS(4,II)*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)+SFKTGRID7(II)
C
C ******************************
C * END OF EMISSION POTENTIALS *
C ******************************
C
      ENDDO
      ENDDO !KKK
      I=I+1 !DAYNUMBER
      ENDDO !L
      CLOSE(21)
      CLOSE(22)

C
C
C ***********************************************
C * INITIALIZATION OF VARIABLES FOR GRAZING ONLY*
C ***********************************************
C
      DO II=1,NXNY
         GRAZING(II)=0
         GRASSGROWTH(II)=0
      DO JJJ=1,7*24
      IGRASSGROWTH(JJJ,II)=0.
      ENDDO
      ENDDO
C

C      
      TIME1=0
C
   21 CONTINUE
C
      WRITE(*,*)'CALCULATION EMISSION VARIATIONS'
      TIME1=0
      OPEN(21,FILE=
     &'/nfs2data03/cas/dyn.pollen.emis/input/W10FIELD.'//CYEAR//'.N2',
     &TYPE='OLD',FORM='UNFORMATTED')
      OPEN(22,FILE=
     &'/nfs2data03/cas/dyn.pollen.emis/input/T2FIELD.'//CYEAR//'.N2',
     &TYPE='OLD',FORM='UNFORMATTED')
      DO M=1,12
      DO KK=1,IMONTH(M)
      DO KKK=1,24
         TIME1=TIME1+1
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      READ(21) U10GRID
      READ(21) V10GRID
      READ(22) TEMPGRID
      DO III=1,NXNY
         VH10GRID(III)=SQRT(U10GRID(III)**2.+V10GRID(III)**2.)
c         VH10GRID(II)=5.
         TEMPGRID(III)=TEMPGRID(III)-273.14
      ENDDO
      DO II=1,NXNY
      IF (TEMPGRID(II).LT.0) THEN
      VENTFACT=0.2
      TSTALD=MAX(0.,(18.0+0.5*(TEMPGRID(I)-0.)))
      ELSE IF (TEMPGRID(II).GE.0.AND.TEMPGRID(II).LT.12.5) THEN
      TSTALD=18.
      VENTFACT=(0.2+(0.38-0.2)*(MAX(0.,TEMPGRID(II))-0.)/(12.5-0.))
      ELSE
      TSTALD=18.+(TEMPGRID(II)-12.5)*0.77
      VENTFACT=0.38 !SAME AS ABOVE
      ENDIF
      FKT1GRID(II)=
     &     (VENTFACT**0.26)*(TSTALD**0.89)*
     &     (1.0*GEMEMIS(1,II)/(SUMISOLATED(II)))
C
C ***************************************
C * EMISSION VARIATION FOR OPEN STABLES *
C ***************************************
C
      IF (TEMPGRID(II).LT.1.) THEN
      TSTALD=4.
      ELSE
      TSTALD=(TEMPGRID(II)+3.)
      ENDIF
      VENTFACT=0.214+
     & 0.014*SIN((243.*TIME1+(REAL(TIME1)))*2.*PI/(24.*REAL(NDAYS)))
       FKT2GRID(II)=
     &      (VENTFACT**0.26)*(TSTALD**0.89)*
     &     (GEMEMIS(2,II)/(SUMNONISOLATED(II)))
C
C **********************************************
C * EMISSION VARIATION FOR STORAGE FASCILITIES *
C **********************************************
C
      IF (TEMPGRID(II).LT.1.) THEN
      TSTALD=1.
      ELSE
      TSTALD=TEMPGRID(II)
      ENDIF
      W10GRID(II)=VH10GRID(II)
      FKT3GRID(II)=(W10GRID(II)**0.26)*(TSTALD**0.89)*
     &     (GEMEMIS(3,II)/(SUMSTORAGE(II)))
C
      W10HELPGRID(II)=EXP(0.0419*W10GRID(II))
C
C ******************************************************
C * EMISSION VARIATION FOR MANURE TYPE 1: EARLY SPRING *
C ******************************************************
C
      C1=(DAYGRIDTEMP(1,II))*24.
      D1=5.*24.
      D1=9.*24.
      WRITE(CMONTH,'(I2.2)') M

      FKT8GRID(II)=EXP(0.0223*TEMPGRID(II))*
     &      W10HELPGRID(II)*
     &     (1./(D1*SQRTTWOPI))*(GEMEMIS(8,II)/(SFKTGRID8(II)))*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)
C
C *****************************************************
C * EMISSION VARIATION FOR MANURE TYPE 2: LATE SPRING *
C *****************************************************
C
      C1=(DAYGRIDTEMP(3,II))*24.
      D1=5.*24.
      D1=9.*24.
      FKT9GRID(II)=EXP(0.0223*TEMPGRID(II))*
     &      W10HELPGRID(II)*
     &     (1./(D1*SQRTTWOPI))*(GEMEMIS(9,II)/(SFKTGRID9(II)))*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)
C
C ************************************************
C * EMISSION VARIATION FOR MANURE TYPE 3: SUMMER *
C ************************************************
C
      C1=(DAYGRIDTEMP(7,II))*24.
      D1=16.*24.
      FKT10GRID(II)=EXP(0.0223*TEMPGRID(II))*
     &      W10HELPGRID(II)*
     &     (1./(D1*SQRTTWOPI))*(GEMEMIS(10,II)/(SFKTGRID10(II)))*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)
C
C
C ************************************************
C * EMISSION VARIATION FOR MANURE TYPE 4: AUTUMN *
C ************************************************
C
      C1=270.*24.
      D1=16.*24.
      FKT11GRID(II)=EXP(0.0223*TEMPGRID(II))*
     &       W10HELPGRID(II)*
     &     (1./(D1*SQRTTWOPI))*(GEMEMIS(11,II)/(SFKTGRID11(II)))*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)
C
C *****************************************************
C * EMISSION VARIATION FOR MANURE TYPE 4a: EMPTY TANKS *
C *****************************************************
C
      C1=270.*24.
      D1=09.*24.
      FKT11aGRID(II)=EXP(0.0223*TEMPGRID(II))*
     &      W10HELPGRID(II)*
     &     (1./(D1*SQRTTWOPI))*(GEMEMIS(12,II)/(SFKTGRID11a(II)))*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)
C
C ****************************************************************************** 
C * EMISSION VARIATION FOR GRAZING CATTLE: ASSUMED TO FOLLOW THE GRASS PATTERN *
C ******************************************************************************
C
      C1=DAYGRIDTEMP(7,II)*24.
      D1=60.*24.
      D1=(GRAZING_GRID(II)/2)*24 !ASSUMING 95% COVERS THE TABULAR VALUE
      FKT14GRID(II)=(1./(D1*SQRTTWOPI))*
     &      EXP(0.0223*TEMPGRID(II))*W10HELPGRID(II)*
     &      (GEMEMIS(14,II)/(SFKTGRID14(II)))*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)
cc      IF (TEMPGRID(II).LT.1.) THEN
cc      TSTALD=1.
cc      ELSE
cc      TSTALD=TEMPGRID(II)
cc      ENDIF
cc      IF (TEMPGRID(II).LT.0.) THEN
cc      GRAZING(II)=0
cc      GRASSGROWTH(II)=0
cc      ELSE
cc      GRASSGROWTH(II)=GRASSGROWTH(II)+TEMPGRID(II)-IGRASSGROWTH(1,II)
cc      ENDIF
C
C ******************************************************************
C * ACCUMULATES GRASSGROWTH LAST 7 DAYS - THE GRAZING REQUIREMENT *
C ******************************************************************
C
cc      DO JJJ=2,7*24
cc      IGRASSGROWTH(JJJ-1,II)=IGRASSGROWTH(JJJ,II)
cc      ENDDO
cc      IGRASSGROWTH(7*24,II)=TEMPGRID(II)
cc
cc      IF (GRASSGROWTH(II).GE.500.AND.
cc     &   (TIME1.GT.15*24.AND.TIME1.LT.350*24)) THEN
cc      GRAZING(II)=1
cc      ELSE
cc      GRAZING(II)=0
cc      ENDIF
cc      FKT14GRID(II)=(GEMEMIS(14,II)/(SFKTGRID14(II)))*
cc     &               GRAZING(II)*(W10GRID(II)**0.26)*TSTALD**0.89

C
C
C ***************************************
C * EMISSION VARIATION FOR WINTER CROPS *
C ***************************************
C
      C1=DAYGRIDTEMP(5,II)*24.
      D1=39.*24.
      FKT4GRID(II)=(1./(D1*SQRTTWOPI))*(GEMEMIS(4,II)/(SFKTGRID4(II)))*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)
C
C
C ***************************************
C * EMISSION VARIATION FOR SPRING CROPS *
C ***************************************
C
      C1=DAYGRIDTEMP(6,II)*24.
      D1=25.*24.
      FKT5GRID(II)=(1./(D1*SQRTTWOPI))*(GEMEMIS(5,II)/(SFKTGRID5(II)))*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)
C
C *********************************************
C * EMISSION VARIATION FOR SPRING SUGAR BEETS *
C *********************************************
C
      C1=DAYGRIDTEMP(8,II)*24.
      D1=45.*24.
      FKT6GRID(II)=(1./(D1*SQRTTWOPI))*(GEMEMIS(6,II)/(SFKTGRID6(II)))*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)
C
C ********************************************
C * EMISSION VARIATION FOR GRASS IN ROTATION *
C ********************************************
C
      C1=DAYGRIDTEMP(7,II)*24.
      D1=60.*24.
      FKT7GRID(II)=(1./(D1*SQRTTWOPI))*(GEMEMIS(7,II)/(SFKTGRID7(II)))*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)
C
C
C *****************************************************
C * EMISSION VARIATION FOR MINERAL FERTILIZER: SPRING *
C *****************************************************
C
      C1=DAYGRIDTEMP(2,II)*24.
      D1=9.*24.
      FKT12GRID(II)=EXP(0.0223*TEMPGRID(II))*
     &      W10HELPGRID(II)* 
     &     (1./(D1*SQRTTWOPI))*(0.9*GEMEMIS(13,II)/(SFKTGRID12(II)))*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)
C
C *****************************************************
C * EMISSION VARIATION FOR MINERAL FERTILIZER: SPRING *
C *****************************************************
C
      C1=DAYGRIDTEMP(7,II)*24.
      D1=16.*24.
      FKT13GRID(II)=EXP(0.0223*TEMPGRID(II))*
     &     W10HELPGRID(II)* 
     &     (1./(D1*SQRTTWOPI))*(0.1*GEMEMIS(13,II)/(SFKTGRID13(II)))*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)
C ************************************************
C * AMMONIA LOSSES FROM NH3 THREATMEANT OF GRASS *
C ************************************************
C
      C1=(DAYGRIDTEMP(4,II)+15)*24.
      D1=9.*24.
      FKT15GRID(II)=GEMEMIS(15,II)*(1./(D1*SQRTTWOPI))*
     &      EXP(-0.5*((TIME1-C1)/D1)**2.)
C
C
C ************************
C * AMMONIA FROM TRAFFIC *
C ************************
C
      FKT16GRID(II)=GEMEMIS(16,II)/(24.*REAL(NDAYS))
C
      ENDDO !NXNY LOOP
      WRITE(*,*) TIME1, GEMEMIS(01,6277), FKT1GRID(6277)
      WRITE(14,'(F6.0,17F20.10)') TIME1, FKT1GRID(6277), FKT2GRID(6277),
     &                                FKT3GRID(6277), FKT4GRID(6277),
     &                                FKT5GRID(6277), FKT6GRID(6277),
     &                                FKT7GRID(6277), FKT8GRID(6277),
     &                                FKT9GRID(6277), FKT10GRID(6277),
     &                                FKT11GRID(6277), FKT11aGRID(6277),
     &                                FKT12GRID(6277), FKT13GRID(6277),
     &                                FKT14GRID(6277), FKT15GRID(6277),
     &                                FKT16GRID(6277)


c      GOTO 22
      IF (DUMB) THEN
      TYPE='FKT01a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT1GRID)
      CLOSE(12)

      TYPE='FKT02a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT2GRID)
      CLOSE(12)
C
      TYPE='FKT03a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT3GRID)
      CLOSE(12)

      TYPE='FKT08a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT8GRID)
      CLOSE(12)
C
      TYPE='FKT09a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'

      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT9GRID)
      CLOSE(12)

      TYPE='FKT10a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT10GRID)
      CLOSE(12)
C
      TYPE='FKT11a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT11GRID)
      CLOSE(12)
C
      TYPE='FKT11b'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT11aGRID)
      CLOSE(12)
C
      TYPE='FKT14a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT14GRID)
      CLOSE(12)
C 
      TYPE='FKT04a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT4GRID)
      CLOSE(12)
C
      TYPE='FKT05a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT5GRID)
      CLOSE(12)

      TYPE='FKT06a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT6GRID)
      CLOSE(12)

      TYPE='FKT07a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT7GRID)
      CLOSE(12)

      TYPE='FKT12a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT12GRID)
      CLOSE(12)

      TYPE='FKT13a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT13GRID)
      CLOSE(12)
C
      TYPE='FKT15a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT15GRID)
      CLOSE(12)

      TYPE='FKT16a'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKT16GRID)
      CLOSE(12)
C DUMBING SUM OF ALL EMISSIONS
      DO II=1,NXNY
         FKTALLGRID(II)=FKT1GRID(II)+
     &   FKT2GRID(II) + FKT3GRID(II) + FKT4GRID(II) +
     &   FKT5GRID(II) + FKT6GRID(II) + FKT7GRID(II) +
     &   FKT8GRID(II) + FKT9GRID(II) + FKT10GRID(II) +
     &   FKT11GRID(II) + FKT11aGRID(II) + FKT12GRID(II) +
     &   FKT13GRID(II) + FKT14GRID(II) + FKT15GRID(II) + 
     &   FKT16GRID(II)
      ENDDO
      TYPE='FKTall'
      WRITE(FILENAME,'(A6,A5,A4,3I2.2,A7)')
     &      TYPE,'.dyn.',
     &      CYEAR,M,KK,KKK,'_d3.bin'
      OPEN(12,FILE=
     &'/nfs2data03/cas/NMR-NH3/NH3emis/'//CYEAR//'/'
     & //CMONTH//'/'//
     &     FILENAME,TYPE='UNKNOWN',FORM='UNFORMATTED')
      CALL PACKDATA(12,FKTALLGRID)
      CLOSE(12)

      ELSE
   22 CONTINUE
      ENDIF

      ENDDO
      ENDDO
      ENDDO
      OPEN(21,FILE='sumgrazing.out')
      DO II=1,NXNY
         WRITE(21,*) II,SUMGRAZING(II)
      ENDDO
      CLOSE(21)

      STOP
C
C
      CLOSE(13)
C
      END
C
C ************************************************************************
C *                       END OF MAIN PROGRAM                            *
C ************************************************************************
C
C
      subroutine packdata(iunit,field)
      integer nx,ny,iunit
      parameter (nx=96,ny=96)
      real field(nx*ny), rmax,rmin,rscale,rmaxint
      integer*2 i2(96*96)
      rmax=field(1)
      rmin=field(1)
      rrmax=0
c$omp parallel do default(shared)
c$omp&private(I) reduction( max:rmax) reduction( min:rmin)
      do i=1,NX*NY
        rmax=max(rmax,field(i))
        rmin=min(rmin,field(i))
      enddo
      rmaxint=32700.
      rscale=(rmax-rmin)/(2*rmaxint)
      if(rscale.eq.0.) rscale=1.
c$omp parallel do default(shared)
c$omp&private(I)
      do i=1,NX*NY
        i2(i)=nint((field(i)-(rmax+rmin)/2)/rscale)
      enddo
      write(iunit) rscale,(rmax+rmin)/2
      write(iunit) i2
      return
      end

