|
SURFEX v7.3
General documentation of Surfex
|
|
SURFEX v7.3
General documentation of Surfex
|
00001 ! ######### 00002 SUBROUTINE CH_EMISSION_FLUX_n(HPROGRAM,PSIMTIME,PSFSV, PRHOA, PTSTEP, KNBTS_MAX) 00003 ! ###################################################################### 00004 !! 00005 !!*** *CH_EMISSION_FLUX_n* - 00006 !! 00007 !! PURPOSE 00008 !! ------- 00009 !! Return a time-dependent emission flux based on tabulated values 00010 !! 00011 !!** METHOD 00012 !! ------ 00013 !! 00014 !! AUTHOR 00015 !! ------ 00016 !! D. Gazen 00017 !! 00018 !! MODIFICATIONS 00019 !! ------------- 00020 !! Original 08/02/00 00021 !! C. Mari 30/10/00 call to MODD_TYPE_EFUTIL and MODD_CST 00022 !! D.Gazen 01/12/03 change emissions handling for surf. externalization 00023 !! P.Tulet 01/01/04 change emission conversion factor 00024 !! P.Tulet 01/01/05 add dust, orilam 00025 !! 00026 !! EXTERNAL 00027 !! -------- 00028 !! 00029 !! 00030 !! IMPLICIT ARGUMENTS 00031 !! ------------------ 00032 USE MODD_SV_n, ONLY: CSV,NSV_CHSBEG,NSV_CHSEND, NSV_AERBEG, NSV_AEREND 00033 USE MODD_TYPE_EFUTIL, ONLY: EMISSVAR_T, PRONOSVAR_T 00034 USE MODD_CSTS, ONLY: NDAYSEC 00035 USE MODD_CH_EMIS_FIELD_n, ONLY: TSEMISS, TSPRONOSLIST, XTIME_SIMUL 00036 USE MODD_CH_SURF_n, ONLY: XCONVERSION 00037 ! 00038 USE MODI_READ_SURF 00039 USE MODI_INIT_IO_SURF_n 00040 USE MODI_END_IO_SURF_n 00041 USE MODI_GET_LUOUT 00042 !UPG*AERO1 00043 USE MODD_CHS_AEROSOL, ONLY: LCH_AERO_FLUX 00044 USE MODI_CH_AER_EMISSION 00045 !UPG*AERO1 00046 !! 00047 !------------------------------------------------------------------------------ 00048 ! 00049 !* 0. DECLARATIONS 00050 ! ----------------- 00051 ! 00052 USE YOMHOOK ,ONLY : LHOOK, DR_HOOK 00053 USE PARKIND1 ,ONLY : JPRB 00054 ! 00055 USE MODI_ABOR1_SFX 00056 ! 00057 IMPLICIT NONE 00058 ! 00059 !* 0.1 declaration of arguments 00060 ! 00061 CHARACTER(LEN=6), INTENT(IN) :: HPROGRAM ! program calling surf. schemes 00062 REAL, INTENT(IN) :: PSIMTIME ! time of simulation in sec UTC 00063 ! (counting from midnight of 00064 ! the current day) 00065 REAL,DIMENSION(:,:), INTENT(INOUT) :: PSFSV ! emission flux in ppp*m/s 00066 REAL, DIMENSION(:), INTENT(IN) :: PRHOA ! air density (kg/m3) 00067 REAL, INTENT(IN) :: PTSTEP ! atmospheric time-step (s) 00068 INTEGER, INTENT(IN) :: KNBTS_MAX !max size of TEMISS%NETIMES 00069 00070 ! 00071 !* 0.2 declaration of local variables 00072 ! 00073 INTEGER :: IVERB ! verbosity level 00074 INTEGER :: KSIZE1D ! 1D size = X*Y physical domain 00075 INTEGER :: JI ! loop control 00076 REAL :: ZALPHA ! interpolation weight 00077 ! 00078 INTEGER :: INBTS ! Number of emission times for a species 00079 INTEGER :: ITIM1,ITIM2 ! first/last time for interpolation 00080 INTEGER :: INDX1,INDX2 ! first/next index for data interpolation 00081 INTEGER :: ISIMTIME, ITPERIOD 00082 CHARACTER (LEN=16) :: YRECFM ! LFI article name 00083 TYPE(PRONOSVAR_T),POINTER :: CURPRONOS !Current pronostic variable 00084 ! 00085 !* 0.3 declaration of saved local variables 00086 ! 00087 CHARACTER(LEN=6), DIMENSION(:), POINTER :: CNAMES 00088 REAL,DIMENSION(SIZE(PSFSV,1),KNBTS_MAX) :: ZWORK ! temporary array for reading data 00089 REAL,DIMENSION(SIZE(PSFSV,1),SIZE(PSFSV,2)) :: ZEMIS ! interpolated in time emission flux 00090 REAL,DIMENSION(SIZE(PSFSV,1)) :: ZFCO ! CO flux 00091 INTEGER :: INEQ ! number of chemical var 00092 !(=NEQ (chimie gaz) + NSV_AER (chimie aerosol) 00093 INTEGER :: IWS ! window size 00094 INTEGER :: IRESP ! return code for I/O 00095 INTEGER :: ILUOUT ! Outputlisting unit 00096 LOGICAL :: LIOINIT ! True if I/O init done 00097 INTEGER :: JW 00098 INTEGER :: ITIME 00099 LOGICAL :: GCO = .FALSE. ! switch if CO emission are available 00100 REAL(KIND=JPRB) :: ZHOOK_HANDLE 00101 ! 00102 !------------------------------------------------------------------------------ 00103 ! 00104 !* EXECUTABLE STATEMENTS 00105 ! --------------------- 00106 ! 00107 IF (LHOOK) CALL DR_HOOK('CH_EMISSION_FLUX_N',0,ZHOOK_HANDLE) 00108 CALL GET_LUOUT(HPROGRAM,ILUOUT) 00109 LIOINIT = .FALSE. 00110 IVERB = 5 00111 KSIZE1D = SIZE(PSFSV,1) 00112 INEQ = SIZE(PSFSV,2) 00113 ! 00114 !------------------------------------------------------------------------------ 00115 ! 00116 !* 3. INTERPOLATE SURFACE FLUXES IN TIME IF NEEDED 00117 ! ------------------------------------------------ 00118 ! 00119 IF (XTIME_SIMUL == 0.) THEN 00120 XTIME_SIMUL = PSIMTIME 00121 ELSE 00122 XTIME_SIMUL = XTIME_SIMUL + PTSTEP 00123 END IF 00124 00125 IF (IVERB >= 5) WRITE(ILUOUT,*) '******** CH_EMISSION_FLUX ********' 00126 DO JI=1,SIZE(TSEMISS) 00127 ! Simulation time (counting from midnight) is saved 00128 ISIMTIME = XTIME_SIMUL 00129 ! 00130 INBTS = SIZE(TSEMISS(JI)%NETIMES) ! 00131 IWS = TSEMISS(JI)%NWS ! Window Size for I/O 00132 INDX1 = TSEMISS(JI)%NDX ! Current data index 00133 ! 00134 IF (INBTS == 1) THEN 00135 ! Time Constant Flux 00136 ! XFWORK already points on data (see ch_buildemiss.f90) 00137 IF (IVERB >= 6) THEN 00138 WRITE(ILUOUT,*) 'NO interpolation for ',TRIM(TSEMISS(JI)%CNAME) 00139 IF (IVERB >= 10 ) WRITE(ILUOUT,*) TSEMISS(JI)%XFWORK 00140 END IF 00141 ELSE 00142 IF (IVERB >= 6) THEN 00143 WRITE(ILUOUT,*) 'Interpolation (T =',ISIMTIME,') : ',TSEMISS(JI)%CNAME 00144 END IF 00145 IF (ISIMTIME < TSEMISS(JI)%NETIMES(1)) THEN 00146 ! Tsim < T(1)=Tmin should not happen but who knows ? 00147 TSEMISS(JI)%NTX = 1 00148 ELSE 00149 ! Check for periodicity when ISIMTIME is beyond last emission time 00150 ! and probably correct ISIMTIME 00151 IF (ISIMTIME > TSEMISS(JI)%NETIMES(INBTS)) THEN 00152 ! Tsim > T(INBTS)=Tmax 00153 ITPERIOD = (1+(TSEMISS(JI)%NETIMES(INBTS)-& 00154 TSEMISS(JI)%NETIMES(TSEMISS(JI)%NPX))/NDAYSEC)*NDAYSEC 00155 ISIMTIME = MODULO(ISIMTIME-TSEMISS(JI)%NETIMES(TSEMISS(JI)%NPX),ITPERIOD)+& 00156 TSEMISS(JI)%NETIMES(TSEMISS(JI)%NPX) 00157 IF (IVERB >= 6) THEN 00158 WRITE(ILUOUT,*) ' ITPERIOD = ', ITPERIOD 00159 WRITE(ILUOUT,*) ' ISIMTIME modifie = ', ISIMTIME 00160 END IF 00161 IF (TSEMISS(JI)%NTX == INBTS .AND. ISIMTIME<TSEMISS(JI)%NETIMES(INBTS)) THEN 00162 ! Update time index NTX 00163 TSEMISS(JI)%NTX = TSEMISS(JI)%NPX 00164 ! Increment data index NDX : NDX correction will occur later 00165 ! to assure 1 <= NDX <= IWS 00166 INDX1 = INDX1 + 1 00167 END IF 00168 END IF 00169 ! 00170 ! search NTX such that : ETIMES(NTX) < ISIMTIME <= ETIMES(NTX+1) 00171 ! and make NDX follow NTX : NDX correction will occur later 00172 ! to assure 1 <= NDX <= IWS 00173 DO WHILE (TSEMISS(JI)%NTX < INBTS) 00174 IF (ISIMTIME >= TSEMISS(JI)%NETIMES(TSEMISS(JI)%NTX+1)) THEN 00175 TSEMISS(JI)%NTX = TSEMISS(JI)%NTX + 1 00176 INDX1 = INDX1 + 1 00177 INDX2 = INDX1 + 1 00178 ELSE 00179 EXIT 00180 END IF 00181 END DO 00182 END IF 00183 ! 00184 ! Check availability of data within memory Window (XEMISDATA(:,1:IWS)) 00185 IF (INDX1 >= IWS) THEN 00186 ! 00187 ! Data index reached the memory window limits 00188 ! 00189 IF (TSEMISS(JI)%LREAD) THEN 00190 ! 00191 ! File must be read to update XEMISDATA array for this species 00192 ! 00193 IF (.NOT. LIOINIT) THEN 00194 ! Must be done once before reading 00195 CALL INIT_IO_SURF_n(HPROGRAM,'FULL ','SURF ','READ ') 00196 IF (IVERB >= 6) WRITE(ILUOUT,*) 'INIT des I/O DONE.' 00197 LIOINIT=.TRUE. 00198 END IF 00199 YRECFM='EMIS_'//TRIM(TSEMISS(JI)%CNAME) 00200 IF (IVERB >= 6)& 00201 WRITE (ILUOUT,*) 'READ emission :',TRIM(YRECFM),& 00202 ', SIZE(ZWORK)=',SIZE(ZWORK,1),INBTS 00203 CALL READ_SURF(HPROGRAM,YRECFM,ZWORK(:,1:INBTS),IRESP) 00204 ! 00205 ! Correction : Replace 999. with 0. value in the Emission FLUX 00206 WHERE(ZWORK(:,1:INBTS) == 999.) 00207 ZWORK(:,1:INBTS) = 0. 00208 END WHERE 00209 WHERE(ZWORK(:,1:INBTS) == 1.E20) 00210 ZWORK(:,1:INBTS) = 0. 00211 END WHERE 00212 DO ITIME=1,INBTS 00213 ZWORK(:,ITIME) = ZWORK(:,ITIME)*XCONVERSION(:) 00214 END DO 00215 ! 00216 ! 00217 IF ((TSEMISS(JI)%NTX+IWS-1) > INBTS) THEN 00218 ! 00219 ! ===== Periodic CASE ===== 00220 ! 00221 IF (IVERB >= 6)& 00222 WRITE (ILUOUT,*) 'Periodic CASE : NPX =',TSEMISS(JI)%NPX 00223 IF (IWS < (INBTS-TSEMISS(JI)%NPX+1)) THEN 00224 ! Window size is smaller then number of periodical times 00225 ! 00226 ! example : IWS=5, NPX=2, INBTS=11, NTX=9 00227 ! NTX NPX 00228 ! | | 00229 ! time index : ...9 10 11 # 2 3 4...11 # 00230 ! old data index :[1 2 3 4 5] 00231 ! new data index : [1 2 3 4 5] 00232 ! | 00233 ! NDX 00234 ! 00235 TSEMISS(JI)%XEMISDATA(:,1:INBTS-TSEMISS(JI)%NTX+1) = & 00236 ZWORK(:,TSEMISS(JI)%NTX:INBTS) 00237 ! 00238 IF (IVERB >= 6) THEN 00239 WRITE(ILUOUT,*) 'Window SIZE smaller than INBTS !' 00240 WRITE(ILUOUT,*) 'Window index, Time index' 00241 DO JW=1,INBTS-TSEMISS(JI)%NTX+1 00242 WRITE(ILUOUT,*) JW,TSEMISS(JI)%NTX+JW-1 00243 END DO 00244 END IF 00245 ! 00246 TSEMISS(JI)%XEMISDATA(:,INBTS-TSEMISS(JI)%NTX+2:IWS) = & 00247 ZWORK(:,TSEMISS(JI)%NPX:TSEMISS(JI)%NPX+IWS-INBTS+TSEMISS(JI)%NTX-2) 00248 ! 00249 IF (IVERB >= 6) THEN 00250 DO JW=INBTS-TSEMISS(JI)%NTX+2,IWS 00251 WRITE(ILUOUT,*) JW,TSEMISS(JI)%NPX+JW-(INBTS-TSEMISS(JI)%NTX+2) 00252 END DO 00253 END IF 00254 INDX1 = 1 00255 INDX2 = 2 00256 ELSE 00257 ! Window size may get smaller AND it will be the last reading 00258 ! 00259 ! example : IWS=6, NPX=7, INBTS=11, NTX=9 00260 ! 00261 ! NTX NPX NTX 00262 ! | | | 00263 ! time index: ...9 10 11 # 7 8 9 10 11 # 00264 ! old data index: ...6] 00265 ! new data index: [1 2 3 4 5] 00266 ! | 00267 ! NDX=NTX-NPX+1 00268 ! 00269 IWS = INBTS-TSEMISS(JI)%NPX+1 00270 TSEMISS(JI)%NWS = IWS 00271 TSEMISS(JI)%XEMISDATA(:,1:IWS) = ZWORK(:,TSEMISS(JI)%NPX:INBTS) 00272 IF (IVERB >= 6) THEN 00273 WRITE(ILUOUT,*) 'Window SIZE equal or greater than INBTS !' 00274 WRITE(ILUOUT,*) 'Window index, Time index' 00275 DO JW=1,IWS 00276 WRITE(ILUOUT,*) JW,TSEMISS(JI)%NPX+JW-1 00277 END DO 00278 END IF 00279 INDX1 = TSEMISS(JI)%NTX-TSEMISS(JI)%NPX+1 00280 INDX2 = MOD((INDX1+1),IWS) 00281 TSEMISS(JI)%LREAD = .FALSE. ! no more reading 00282 END IF 00283 ELSE 00284 ! 00285 ! ===== NON periodic (normal) CASE ===== 00286 ! 00287 ! example : with IWS=5, the window moves forward 00288 ! NTX 00289 ! | 00290 ! time index : 1 2 3 4 5 6 7 8 9 10 11 ... INBTS # 00291 ! old data index :[1 2 3 4 5] 00292 ! new data index : [1 2 3 4 5] 00293 ! | 00294 ! NDX 00295 ! 00296 TSEMISS(JI)%XEMISDATA(:,1:IWS) = ZWORK(:,TSEMISS(JI)%NTX:TSEMISS(JI)%NTX+IWS-1) 00297 IF (IVERB >= 6) THEN 00298 WRITE(ILUOUT,*) 'Window index, Time index' 00299 DO JW=1,IWS 00300 WRITE(ILUOUT,*) JW,TSEMISS(JI)%NTX+JW-1 00301 END DO 00302 END IF 00303 INDX1 = 1 00304 INDX2 = 2 00305 END IF 00306 ELSE 00307 ! Data is already in memory because window size is sufficient 00308 ! to hold INBTS emission times => simply update NDX according to NTX 00309 ! 00310 IF (IWS==INBTS) THEN 00311 ! 00312 ! 'Window size' = 'Nb emis times' at INIT (ch_init_emission) 00313 ! so NDX must be set equal to NTX (the window does not move) 00314 ! example : 00315 ! NPX NTX 00316 ! | | 00317 ! time index : 1 2 3 ... INBTS 00318 ! data index : [1 2 3 ... INBTS] 00319 ! | 00320 ! NDX 00321 00322 INDX1 = TSEMISS(JI)%NTX 00323 INDX2 = INDX1+1 00324 IF (INDX2 > IWS) INDX2=TSEMISS(JI)%NPX 00325 ELSE 00326 ! 00327 ! Windows size changed during periodic case 00328 ! NDX must be equal to NTX - NPX + 1 00329 ! (the window does not move) 00330 ! example : 00331 ! NTX 00332 ! | 00333 ! time index : NPX NPX+1 NPX+2 ... INBTS 00334 ! data index : [1 2 3 ... IWS] 00335 ! | 00336 ! NDX 00337 INDX1 = TSEMISS(JI)%NTX-TSEMISS(JI)%NPX+1 00338 INDX2 = MOD((INDX1+1),IWS) 00339 END IF 00340 END IF 00341 ELSE ! (INDX1 < IWS) 00342 INDX2 = INDX1+1 00343 END IF 00344 ! 00345 ! Don't forget to update NDX with new value INDX1 00346 TSEMISS(JI)%NDX = INDX1 00347 ! 00348 ! Compute both times for interpolation 00349 IF (TSEMISS(JI)%NTX < INBTS) THEN 00350 ITIM1 = TSEMISS(JI)%NETIMES(TSEMISS(JI)%NTX) 00351 ITIM2 = TSEMISS(JI)%NETIMES(TSEMISS(JI)%NTX+1) 00352 ELSE 00353 ITIM1 = TSEMISS(JI)%NETIMES(INBTS) 00354 ITIM2 = TSEMISS(JI)%NETIMES(TSEMISS(JI)%NPX)+ITPERIOD 00355 END IF 00356 ! 00357 ! Interpolate variables in time -> update XFWORK 00358 ! 00359 ! 00360 ! time : ITIM1...Tsim...ITIM2 00361 ! | | 00362 ! data index : INDX1 INDX2 00363 ! 00364 ! 00365 ZALPHA = (REAL(ISIMTIME) - ITIM1) / (ITIM2-ITIM1) 00366 TSEMISS(JI)%XFWORK(:) = ZALPHA*TSEMISS(JI)%XEMISDATA(:,INDX2) +& 00367 (1.-ZALPHA)*TSEMISS(JI)%XEMISDATA(:,INDX1) 00368 IF (IVERB >= 6) THEN 00369 WRITE(ILUOUT,*) ' Current time INDEX : ',TSEMISS(JI)%NTX 00370 WRITE(ILUOUT,*) ' TIME : ',ISIMTIME, ' (',ITIM1,',',ITIM2,')' 00371 WRITE(ILUOUT,*) ' Window size : ',TSEMISS(JI)%NWS 00372 WRITE(ILUOUT,*) ' Current data INDEX : ',INDX1,INDX2 00373 IF (IVERB >= 10) WRITE(ILUOUT,*) ' FLUX : ',TSEMISS(JI)%XFWORK 00374 END IF 00375 END IF 00376 END DO 00377 ! 00378 ! Agregation : flux computation 00379 ! 00380 ZEMIS(:,:) = 0. 00381 ! 00382 ! Point on head of Pronostic variable list 00383 ! to cover the entire list. 00384 IF (NSV_AEREND > 0) THEN 00385 CNAMES=>CSV(NSV_CHSBEG:NSV_AEREND) 00386 ELSE 00387 CNAMES=>CSV(NSV_CHSBEG:NSV_CHSEND) 00388 END IF 00389 CURPRONOS=>TSPRONOSLIST 00390 DO WHILE(ASSOCIATED(CURPRONOS)) 00391 IF (CURPRONOS%NAMINDEX > INEQ) THEN 00392 WRITE(ILUOUT,*) 'FATAL ERROR in CH_EMISSION_FLUXN : SIZE(ZEMIS,2) =',& 00393 INEQ,', INDEX bugge =',CURPRONOS%NAMINDEX 00394 CALL ABOR1_SFX('CH_EMISSION_FLUXN: FATAL ERROR') 00395 END IF 00396 00397 ZEMIS(:,CURPRONOS%NAMINDEX) = 0. 00398 ! 00399 ! Loop on the number of agreg. coeff. 00400 DO JI=1,CURPRONOS%NBCOEFF 00401 ! Compute agregated flux 00402 ZEMIS(:,CURPRONOS%NAMINDEX) = ZEMIS(:,CURPRONOS%NAMINDEX)+& 00403 CURPRONOS%XCOEFF(JI)*TSEMISS(CURPRONOS%NEFINDEX(JI))%XFWORK(:) 00404 END DO 00405 00406 IF (IVERB >= 6) THEN 00407 WRITE(ILUOUT,*) 'Agregation for ',CNAMES(CURPRONOS%NAMINDEX) 00408 IF (IVERB >= 10) WRITE(ILUOUT,*) 'ZEMIS = ',ZEMIS(:,CURPRONOS%NAMINDEX) 00409 END IF 00410 IF ((CNAMES(CURPRONOS%NAMINDEX) == "CO") .AND. ANY(ZEMIS(:,CURPRONOS%NAMINDEX).GT.0.)) THEN 00411 ZFCO(:) = ZEMIS(:,CURPRONOS%NAMINDEX) 00412 GCO = .TRUE. 00413 END IF 00414 00415 CURPRONOS=>CURPRONOS%NEXT 00416 ! 00417 END DO 00418 ! 00419 IF ((LCH_AERO_FLUX).AND.(NSV_AERBEG > 0)) THEN 00420 IF (GCO) THEN 00421 CALL CH_AER_EMISSION(ZEMIS, PRHOA, CSV, NSV_CHSBEG, PFCO=ZFCO) 00422 ELSE 00423 CALL CH_AER_EMISSION(ZEMIS, PRHOA, CSV, NSV_CHSBEG) 00424 ENDIF 00425 END IF 00426 ! 00427 PSFSV(:,:) = PSFSV(:,:) + ZEMIS(:,:) 00428 ! 00429 IF (LIOINIT) CALL END_IO_SURF_n(HPROGRAM) 00430 ! 00431 IF (IVERB >= 6) WRITE(ILUOUT,*) '******** END CH_EMISSION_FLUX ********' 00432 IF (LHOOK) CALL DR_HOOK('CH_EMISSION_FLUX_N',1,ZHOOK_HANDLE) 00433 ! 00434 END SUBROUTINE CH_EMISSION_FLUX_n
1.8.0
