!MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt !MNH_LIC for details. version 1. !------------------------------------------------------------------------------ ! ##################### MODULE MODI_EMIS_INIT ! ##################### ! INTERFACE ! SUBROUTINE EMIS_INIT(KNUMB) ! INTEGER ,INTENT(IN) :: KNUMB ! END SUBROUTINE EMIS_INIT ! END INTERFACE ! END MODULE MODI_EMIS_INIT ! ! ########################### SUBROUTINE EMIS_INIT(KNUMB) ! ########################### ! !! !!**** EMIS_INIT stores the pollutant emissions !! !! PURPOSE !! ------- !!**** The purpose of this routine is to read an ASCII file containing !! the information of the pollutant emissions and to store the information !! in a dedicated array !! the storage is split in space (parallel computation) and in time (user choice) !! !! METHOD !! ------ !!**** Three main steps !! - read input ASCII files !! - eliminate the information used out the processor domain !! - class the retained information in time !! units : x(m),y(m),z(m),t(s),drhodt(kg.m-3.s-1) !! !! EXTERNAL !! -------- !! SUBROUTINE ? !! !! IMPLICIT ARGUMENTS !! ------------------ !! MODD_? !! !! REFERENCE !! --------- !! !! AUTHOR !! ------ !! Franck Auguste (CERFACS-AE) !! !! MODIFICATIONS !! ------------- !! Original 01/10/2017 !! !------------------------------------------------------------------------------ ! !**** 0. DECLARATIONS ! --------------- ! ! module USE MODE_POS USE MODE_ll USE MODE_IO_ll USE MODD_VAR_ll, ONLY: IP USE MODD_ARGSLIST_ll, ONLY : LIST_ll USE MODD_TIME_n, ONLY: TDTCUR ! ! declaration USE MODD_EMIS_PARAM_n USE MODD_FIELD_n USE MODD_DIM_n, ONLY: NIMAX,NJMAX,NKMAX USE MODD_PARAMETERS, ONLY: JPVEXT,JPHEXT USE MODD_GRID USE MODD_GRID_n, ONLY: XXHAT,XYHAT,XZZ USE MODD_METRICS_n, ONLY: XDXX,XDYY,XDZZ,XDZX,XDZY USE MODD_CONF, ONLY: NHALO,CCONF USE MODD_SUB_MODEL_n USE MODD_REF_n ! ! interface USE MODD_GRID USE MODD_CST USE MODD_GRID_n USE MODE_GRIDPROJ USE MODI_SECOND_MNH USE MODN_PARAM_n ! IMPLICIT NONE ! !------------------------------------------------------------------------------ ! ! 0.1 declarations of arguments ! INTEGER ,INTENT(IN) :: KNUMB ! number of input files ! !------------------------------------------------------------------------------ ! ! 0.2 declaration of local variables ! INTEGER :: IIB,IIE,IJB,IJE,IKB,IKE ! proc size INTEGER :: JN,JM ! loop index INTEGER :: ILUEMI1,IRESPEMI1,ILUEMI2,IRESPEMI2 ! integers for open/read files INTEGER :: INUM11,INUM12 ! line number and species number REAL :: ZDELTX1,ZDELTT1 ! space and time step REAL :: ZLATREF1,ZLONREF1,ZXREF1,ZYREF1,ZZREF1,ZZREF0 ! reference location REAL :: ZTREF1,ZTREF0 ! time location REAL :: ZXMIN,ZYMIN,ZXMAX,ZYMAX,ZZMIN,ZZMAX ! proc location REAL, DIMENSION(:,:), ALLOCATABLE :: PTMP1 ! first writing of the input file INTEGER, DIMENSION(:,:), ALLOCATABLE :: INUMB1A ! searching the array size REAL :: ZXTMP0,ZYTMP0,ZZTMP0,ZTTMP0 ! temporary time/space location INTEGER :: ITEMP1,ITEMP2,ITEMP3 ! temporary integer INTEGER :: IMAXI1A ! the array size REAL :: ZTIME0,ZTIME1, ZTIME2 ! estimate CPU time REAL :: ZCOEFZ ! ground weighting INTEGER :: ILUOUT INTEGER :: IINFO_ll REAL :: ZTMP_MIN,ZTMP_MAX,ZTMP_BIS,ZRESTA !------------------------------------------------------------------------------ ! ! XEMIS_DATA(:,:,:,:) ! the output file(24_hour,60_min,nb_data,nb_species) ! NEMIS_NUMB=INUM2 ! nb_species ! NEMIS_MAXI=IMAXIA ! nb_data ! XEMIS_HEIG=? ! reference altitude ! !------------------------------------------------------------------------------ ! !* *** 1. PRELIMINARIES ! ----------------- ! ZRESTA=0. IF (CCONF=='RESTA') ZRESTA=TDTCUR%TIME CALL GET_INDICE_ll(IIB,IJB,IIE,IJE) IKB = 1+JPVEXT IKE = SIZE(XZZ,3)-JPVEXT XT_EMIS_INIT = 0. IF (IP==1) THEN IF (KNUMB==0) THEN WRITE(*,*) '************************' WRITE(*,*) '***** NEMIS = TRUE *****' WRITE(*,*) '*** FILE IS REQUIRED ***' WRITE(*,*) '*** (stop execution) ***' WRITE(*,*) '************************' STOP ENDIF IF (NEMIS_HOUR==0) THEN WRITE(*,*) '************************' WRITE(*,*) '***** NEMIS = TRUE *****' WRITE(*,*) '*** HOUR IS REQUIRED ***' WRITE(*,*) '*** (stop execution) ***' WRITE(*,*) '************************' STOP ENDIF ENDIF ! !* *** 2. READ AND SPLIT ! ------------------ ! CALL SECOND_MNH(ZTIME1) IF (KNUMB>=1) THEN CALL OPEN_ll(UNIT=ILUEMI1 , FILE='emis1.nam', IOSTAT=IRESPEMI1 , FORM='FORMATTED' , & STATUS='OLD ',ACCESS='SEQUENTIAL ', ACTION='READ', MODE = GLOBAL ) ! CALL OPEN_ll(UNIT=ILUEMI2 , FILE='emis2.nam', IOSTAT=IRESPEMI2 , FORM='FORMATTED' , & ! STATUS='OLD ',ACCESS='SEQUENTIAL ', ACTION='READ', MODE = GLOBAL ) READ(UNIT=ILUEMI1,FMT=*) INUM11,INUM12 ! line number and species number (column except x,y,z,t) READ(UNIT=ILUEMI1,FMT=*) ZDELTX1,ZDELTT1 ! space and time step READ(UNIT=ILUEMI1,FMT=*) ZLATREF1,ZLONREF1 ! reference location READ(UNIT=ILUEMI1,FMT=*) ZTREF1,ZZREF1,ZTREF0,ZZREF0! time/altitude location (IESTA:12068s(3.5h) to 81827s(22.5h) and 300m) ! time start simulation + altitude start/restart (max alt emission : 44m) IF (IP==1) THEN IF (INUM12>11) THEN WRITE(*,*) '************************' WRITE(*,*) '***** NEMIS = TRUE *****' WRITE(*,*) '*** NUMB<12 REQUIRED ***' WRITE(*,*) '*** (stop execution) ***' WRITE(*,*) '************************' STOP ENDIF ENDIF CALL SM_XYHAT_S(XLATORI,XLONORI,ZLATREF1,ZLONREF1,ZXREF1,ZYREF1) ! to activate if use of lat/lon coordinates ! ZXREF1=ZLATREF1 ! ZYREF1=ZLONREF1 XEMIS_HEIG = ZZREF1 XEMIS1_DELX = ZDELTX1 XEMIS1_DELT = ZDELTT1 ZXMIN = XXHAT(IIB ) -(1.+XEMIS_RESI)*ZDELTX1 ZYMIN = XYHAT(IJB ) -(1.+XEMIS_RESI)*ZDELTX1 ZZMIN = XZZ(IIB,IJB,IKB)-(0.+XEMIS_RESI)*ZDELTX1 ZXMAX = XXHAT(IIE+1) +(1.+XEMIS_RESI)*ZDELTX1 ZYMAX = XYHAT(IJE+1) +(1.+XEMIS_RESI)*ZDELTX1 ZZMAX = XZZ(IIB,IJE,IKE)+(0.+XEMIS_RESI)*ZDELTX1 ZTMP_BIS=0. ZTMP_MAX=0. ZTMP_MIN=1000000. ! one to eleven variables can be specified ALLOCATE(PTMP1(INUM11,INUM12+4)) PTMP1(:,:) = 0. JM = 0 DO JN=1,INUM11 JM = JM + 1 IF (INUM12== 1) READ(UNIT=ILUEMI1,FMT=*) PTMP1(JN, 1),PTMP1(JN, 2),PTMP1(JN, 3),PTMP1(JN, 4),PTMP1(JN, 5) IF (INUM12== 2) READ(UNIT=ILUEMI1,FMT=*) PTMP1(JN, 1),PTMP1(JN, 2),PTMP1(JN, 3),PTMP1(JN, 4),PTMP1(JN, 5),& PTMP1(JN, 6)!,PTMP1(JN, 7),PTMP1(JN, 8),PTMP1(JN, 9),PTMP1(JN,10) IF (INUM12== 3) READ(UNIT=ILUEMI1,FMT=*) PTMP1(JN, 1),PTMP1(JN, 2),PTMP1(JN, 3),PTMP1(JN, 4),PTMP1(JN, 5),& PTMP1(JN, 6),PTMP1(JN, 7)!,PTMP1(JN, 8),PTMP1(JN, 9),PTMP1(JN,10) IF (INUM12== 4) READ(UNIT=ILUEMI1,FMT=*) PTMP1(JN, 1),PTMP1(JN, 2),PTMP1(JN, 3),PTMP1(JN, 4),PTMP1(JN, 5),& PTMP1(JN, 6),PTMP1(JN, 7),PTMP1(JN, 8)!,PTMP1(JN, 9),PTMP1(JN,10) IF (INUM12== 5) READ(UNIT=ILUEMI1,FMT=*) PTMP1(JN, 1),PTMP1(JN, 2),PTMP1(JN, 3),PTMP1(JN, 4),PTMP1(JN, 5),& PTMP1(JN, 6),PTMP1(JN, 7),PTMP1(JN, 8),PTMP1(JN, 9)!,PTMP1(JN,10) IF (INUM12== 6) READ(UNIT=ILUEMI1,FMT=*) PTMP1(JN, 1),PTMP1(JN, 2),PTMP1(JN, 3),PTMP1(JN, 4),PTMP1(JN, 5),& PTMP1(JN, 6),PTMP1(JN, 7),PTMP1(JN, 8),PTMP1(JN, 9),PTMP1(JN,10) IF (INUM12== 7) READ(UNIT=ILUEMI1,FMT=*) PTMP1(JN, 1),PTMP1(JN, 2),PTMP1(JN, 3),PTMP1(JN, 4),PTMP1(JN, 5),& PTMP1(JN, 6),PTMP1(JN, 7),PTMP1(JN, 8),PTMP1(JN, 9),PTMP1(JN,10),& PTMP1(JN,11)!,PTMP1(JN,12),PTMP1(JN,13),PTMP1(JN,14),PTMP1(JN,15) IF (INUM12== 8) READ(UNIT=ILUEMI1,FMT=*) PTMP1(JN, 1),PTMP1(JN, 2),PTMP1(JN, 3),PTMP1(JN, 4),PTMP1(JN, 5),& PTMP1(JN, 6),PTMP1(JN, 7),PTMP1(JN, 8),PTMP1(JN, 9),PTMP1(JN,10),& PTMP1(JN,11),PTMP1(JN,12)!,PTMP1(JN,13),PTMP1(JN,14),PTMP1(JN,15) IF (INUM12== 9) READ(UNIT=ILUEMI1,FMT=*) PTMP1(JN, 1),PTMP1(JN, 2),PTMP1(JN, 3),PTMP1(JN, 4),PTMP1(JN, 5),& PTMP1(JN, 6),PTMP1(JN, 7),PTMP1(JN, 8),PTMP1(JN, 9),PTMP1(JN,10),& PTMP1(JN,11),PTMP1(JN,12),PTMP1(JN,13)!,PTMP1(JN,14),PTMP1(JN,15) IF (INUM12==10) READ(UNIT=ILUEMI1,FMT=*) PTMP1(JN, 1),PTMP1(JN, 2),PTMP1(JN, 3),PTMP1(JN, 4),PTMP1(JN, 5),& PTMP1(JN, 6),PTMP1(JN, 7),PTMP1(JN, 8),PTMP1(JN, 9),PTMP1(JN,10),& PTMP1(JN,11),PTMP1(JN,12),PTMP1(JN,13),PTMP1(JN,14)!,PTMP1(JN,15) IF (INUM12==11) READ(UNIT=ILUEMI1,FMT=*) PTMP1(JN, 1),PTMP1(JN, 2),PTMP1(JN, 3),PTMP1(JN, 4),PTMP1(JN, 5),& PTMP1(JN, 6),PTMP1(JN, 7),PTMP1(JN, 8),PTMP1(JN, 9),PTMP1(JN,10),& PTMP1(JN,11),PTMP1(JN,12),PTMP1(JN,13),PTMP1(JN,14),PTMP1(JN,15) ! ZTMP_MIN = MIN(PTMP1(JN, 4),ZTMP_MIN) ! ZTMP_MAX = MAX(PTMP1(JN, 4),ZTMP_MAX) ! IF ((PTMP1(JN, 1)>0.).AND.(PTMP1(JN, 2)>0.).AND.(PTMP1(JN, 1)<6000.).AND.(PTMP1(JN, 2)<4000.)) ZTMP_BIS = MAX(PTMP1(JN, 3),ZTMP_BIS) ! IF (MOD(JM,10)==0.AND.(IP==1)) WRITE(*,*)JM,ZTMP_BIS ! IF (MOD(JM,1000)==0.AND.(IP==1)) WRITE(*,*)JM,PTMP1(JN,10),PTMP1(JN,11) ! IF (((PTMP1(JN, 5)>5.).OR.(PTMP1(JN, 5)<-5.)).AND.(IP==1)) WRITE(*,*)'=05=',JM ! IF (((PTMP1(JN, 6)>5.).OR.(PTMP1(JN, 6)<-5.)).AND.(IP==1)) WRITE(*,*)'=06=',JM ! IF (((PTMP1(JN, 7)>5.).OR.(PTMP1(JN, 7)<-5.)).AND.(IP==1)) WRITE(*,*)'=07=',JM ! IF (((PTMP1(JN, 8)>5.).OR.(PTMP1(JN, 8)<-5.)).AND.(IP==1)) WRITE(*,*)'=08=',JM ! IF (((PTMP1(JN, 9)>5.).OR.(PTMP1(JN, 9)<-5.)).AND.(IP==1)) WRITE(*,*)'=09=',JM ! IF (((PTMP1(JN,10)>5.).OR.(PTMP1(JN,10)<-5.)).AND.(IP==1)) WRITE(*,*)'=10=',JM ! IF (((PTMP1(JN,11)>5.).OR.(PTMP1(JN,11)<-5.)).AND.(IP==1)) WRITE(*,*)'=11=',JM ENDDO ! searching of the array size after time/space splitting ALLOCATE(INUMB1A(NEMIS_HOUR,60)) INUMB1A=0 IMAXI1A=0 DO JN=1,INUM11 ZXTMP0=ZXREF1+PTMP1(JN, 1)*ZDELTX1 ZYTMP0=ZYREF1+PTMP1(JN, 2)*ZDELTX1 ZZTMP0=ZZREF1+(PTMP1(JN, 3)-ZZREF0)*ZDELTX1-ZZREF0 ZTTMP0=ZTREF1+(PTMP1(JN, 4)-ZTREF0)*ZDELTT1-ZTREF0-ZRESTA ITEMP1=INT((ZTTMP0/3600.)) ITEMP2=INT((ZTTMP0-ITEMP1*3600.)/60.) ITEMP1=1+ITEMP1 ITEMP2=1+ITEMP2 IF (ITEMP1<=0) GOTO 101 IF (ITEMP1>NEMIS_HOUR) GOTO 102 ITEMP2=MAX(1,ITEMP2) ITEMP2=MIN(60,ITEMP2) IF ((ZXTMP0.GT. ZXMIN).AND.(ZXTMP0.LT.ZXMAX).AND. & (ZYTMP0.GT. ZYMIN).AND.(ZYTMP0.LT.ZYMAX).AND. & (ZZTMP0.GT. ZZMIN).AND.(ZZTMP0.LT.ZZMAX)) THEN INUMB1A(ITEMP1,ITEMP2)=INUMB1A(ITEMP1,ITEMP2)+1 ENDIF 101 CONTINUE 102 CONTINUE ENDDO DO JN=1,NEMIS_HOUR DO JM=1,60 IMAXI1A=MAX(IMAXI1A,INUMB1A(JN,JM)) ENDDO ENDDO ! storing the array after time/space splitting ALLOCATE(XEMIS1_DATA(NEMIS_HOUR,60,IMAXI1A,INUM12+4)) XEMIS1_DATA(:,:,:,:)=0. XEMIS1_DATA(:,:,:,1)=-9999. XEMIS1_DATA(:,:,:,2)=-9999. XEMIS1_DATA(:,:,:,3)=-9999. XEMIS1_DATA(:,:,:,4)=-9999. INUMB1A(:,:)=1 DO JN=1,INUM11 ZXTMP0=ZXREF1+PTMP1(JN, 1)*ZDELTX1 ZYTMP0=ZYREF1+PTMP1(JN, 2)*ZDELTX1 ZZTMP0=ZZREF1+(PTMP1(JN, 3)-ZZREF0)*ZDELTX1-ZZREF0 ZTTMP0=ZTREF1+(PTMP1(JN, 4)-ZTREF0)*ZDELTT1-ZTREF0-ZRESTA ZCOEFZ=1. IF (ZZTMP0<-ZZMIN) ZCOEFZ=2. ITEMP1=INT((ZTTMP0/3600.)) ITEMP2=INT((ZTTMP0-ITEMP1*3600.)/60.) ITEMP1=1+ITEMP1 ITEMP2=1+ITEMP2 IF (ITEMP1<=0) GOTO 201 IF (ITEMP1>NEMIS_HOUR) GOTO 202 ITEMP2=MAX(1,ITEMP2) ITEMP2=MIN(60,ITEMP2) IF ((ZXTMP0.GT. ZXMIN).AND.(ZXTMP0.LT.ZXMAX).AND. & (ZYTMP0.GT. ZYMIN).AND.(ZYTMP0.LT.ZYMAX).AND. & (ZZTMP0.GT. ZZMIN).AND.(ZZTMP0.LT.ZZMAX)) THEN ITEMP3 = INUMB1A(ITEMP1,ITEMP2) XEMIS1_DATA(ITEMP1,ITEMP2,ITEMP3,1)=ZXTMP0 XEMIS1_DATA(ITEMP1,ITEMP2,ITEMP3,2)=ZYTMP0 XEMIS1_DATA(ITEMP1,ITEMP2,ITEMP3,3)=ZZTMP0 XEMIS1_DATA(ITEMP1,ITEMP2,ITEMP3,4)=ZTTMP0 XEMIS1_DATA(ITEMP1,ITEMP2,ITEMP3,5:(5+INUM12-1))=ZCOEFZ*PTMP1(JN,5:(5+INUM12-1))/ZDELTX1**3./ZDELTT1 !(kg.m-3.s-1) INUMB1A(ITEMP1,ITEMP2)=INUMB1A(ITEMP1,ITEMP2)+1 ENDIF 201 CONTINUE 202 CONTINUE ENDDO DEALLOCATE(PTMP1) ENDIF ! CALL SECOND_MNH(ZTIME2) XT_EMIS_INIT=ZTIME2-ZTIME1 ! RETURN ! END SUBROUTINE EMIS_INIT