! ! $Id: stratemit.F90 2022-07-04 mmarchand $ ! SUBROUTINE STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,tr_seri,& m_emiss_vol_daily,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_spec,& stretchlong,ispeci,id_species_total) USE dimphy, ONLY : klon,klev USE strataer_local_var_mod USE phys_cal_mod USE phys_local_var_mod, ONLY: d_q_emiss USE yomcst_mod_h, only : RD, RPI, RG USE geometry_mod, ONLY : cell_area, boundslat USE aerophys USE infotrac_phy USE mod_phys_lmdz_mpi_data, ONLY : is_mpi_root USE mod_grid_phy_lmdz, ONLY: nbp_lon IMPLICIT NONE ! Input argument !--------------- REAL,INTENT(IN) :: pdtphys, pdt ! Pas d'integration pour la physique(seconde) REAL,DIMENSION(klon),INTENT(IN) :: xlat ! latitudes pour chaque point REAL,DIMENSION(klon),INTENT(IN) :: xlon ! longitudes pour chaque point REAL,DIMENSION(klon,klev),INTENT(IN) :: pplay ! pression pour le mileu de chaque couche (en Pa) REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs ! pression pour chaque inter-couche (en Pa) REAL,DIMENSION(klon,klev),INTENT(IN) :: t_seri ! Temperature REAL,INTENT(IN) :: latmin,latmax,lonmin,lonmax REAL,INTENT(IN) :: m_emiss_vol_daily ! daily injection mass emission REAL,INTENT(IN) :: stretchlong ! length to stretch emission (day) REAL,DIMENSION(klon,klev) :: m_air_gridbox ! mass ofair in every grid box[kg REAL,INTENT(IN) :: sigma_alt, altemiss INTEGER,INTENT(IN) :: id_spec, ispeci, id_species_total ! Output argument !---------------- REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT) :: tr_seri ! Concentration Traceur [U/KgA] ! Local variables !---------------- INTEGER :: ieru, i, k, i_int, ilon, ilev REAL :: dlat_loc REAL,DIMENSION(klev+1) :: altLMDz ! altitude of layer interfaces in m REAL :: zrho ! Density of air [kg/m3] REAL,DIMENSION(klev) :: zdm ! mass of atm. model layer in kg REAL :: zdz ! thickness of atm. model layer in m REAL :: f_lay_sum ! sum of layer emission fractions REAL,DIMENSION(klev) :: f_lay_emiss ! emiss fraction for every emiss for every vertical layer REAL :: emission ! emission REAL :: theta_min, theta_max, theta ! for SAI computation between two latitudes IF (is_mpi_root) THEN WRITE(*,*) 'IN STRATEMIT: date from phys_cal_mod=',year_cur,'-',& & mth_cur,'-',day_cur,'-',hour,' flh2o=',flh2o ENDIF !--calculate mass of air in every grid box !NL: DONE ABOVE DIRECTLY IN MAIN LOOP I,K ! DO ilon=1, klon ! DO ilev=1, klev ! m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon) ! ENDDO ! ENDDO DO i=1,klon dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3)) dlon=360./2./FLOAT(nbp_lon) ! dlat = half difference of boundary latitudes IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT, dlat=',dlat_loc,xlat(i),xlon(i) theta_min = max(xlat(i)-dlat_loc,latmin) theta_max = min(xlat(i)+dlat_loc,latmax) IF ( xlat(i).GE.latmin-dlat_loc .AND. & & xlat(i).LT.latmax+dlat_loc .AND. & & xlon(i).GE.lonmin-dlon .AND. & & xlon(i).LT.lonmax+dlon ) THEN ! WRITE(*,*) 'coordinates of volcanic injection point=',& & xlat(i),xlon(i),day_cur,mth_cur,year_cur WRITE(*,*) 'DD m_emiss_vol_daily=', & & m_emiss_vol_daily !compute altLMDz altLMDz(:)=0.0 DO k=1, klev zrho=pplay(i,k)/t_seri(i,k)/RD !air density in kg/m3 zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG !mass of layer in kg zdz=zdm(k)/zrho !thickness of layer in m altLMDz(k+1)=altLMDz(k)+zdz !altitude of interface ENDDO CALL STRATDISTRIB(altLMDz,altemiss,sigma_alt,f_lay_emiss) IF (flag_emit.EQ.3) then theta=(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ & & (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI)) ELSE theta=1. ENDIF !vertically distributed emission DO k=1, klev !--calculate mass of air in every grid box m_air_gridbox(i,k)=(paprs(i,k)-paprs(i,k+1))/RG*cell_area(i) ! stretch emission over stretchlong period emission=m_emiss_vol_daily/m_air_gridbox(i,k)*f_lay_emiss(k)/stretchlong/ & & (86400.-pdt)*theta IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: emission avant/apres', & & 'alt= ', altLMDz(k), & & 'flh2o= ',flh2o, & & 'id_speac= ',id_spec,id_species_total, & & 'emission= ',emission, & & 'pdtphys= ',pdtphys, & & 'rapport m_emiss/m_air*f_lay= ', m_emiss_vol_daily/m_air_gridbox(i,k)*f_lay_emiss(k), & & 'stretchlong= ', stretchlong, & & 'theta= ', theta IF(emission < 1.E-34) emission = 0.0 IF (flh2o.EQ.0) THEN IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: tr_ser avant/apres',& & 'i= ',i,'k= ',k, 'flh2o= ',flh2o, & & tr_seri(i,k,id_spec), & & tr_seri(i,k,id_spec)+emission*pdtphys tr_seri(i,k,id_spec)=tr_seri(i,k,id_spec)+emission*pdtphys IF (id_species_total.NE.0) THEN tr_seri(i,k,id_species_total)=tr_seri(i,k,id_species_total)+emission*pdtphys ENDIF ELSE IF(flh2o.EQ.1) THEN d_q_emiss(i,k)=emission*pdtphys IF(d_q_emiss(i,k) > 1.E34) THEN WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) > 1.E34 = ',d_q_emiss(i,k) ELSE IF(d_q_emiss(i,k) < -1.E34) THEN WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) < -1.E34 = ',d_q_emiss(i,k) ENDIF IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: ',& & 'i= ',i,'k= ',k, 'flh2o= ',flh2o, & & 'emission= ',emission, & & 'd_q_emiss(i,k)= ',d_q_emiss(i,k) IF(d_q_emiss(i,k) > 1.E34) THEN d_q_emiss(i,k) = 1.E34 WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) > 1.E34 : ',d_q_emiss(i,k) ELSE IF(d_q_emiss(i,k) < -1.E34) THEN d_q_emiss(i,k) = -1.E34 WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) < -1.E34 : ',d_q_emiss(i,k) ENDIF ENDIF budg_emi(i,ispeci)=budg_emi(i,ispeci)+emission*zdm(k) ENDDO ENDIF ! emission grid cell ENDDO ! klon loop END SUBROUTINE STRATEMIT