source: LMDZ6/trunk/libf/phylmd/StratAer/stratemit.F90 @ 4601

Last change on this file since 4601 was 4601, checked in by dcugnet, 11 months ago

StratAer? commit, N. Lebas
File size: 7.0 KB
Line 
1!
2! $Id: stratemit.F90  2022-07-04 mmarchand $
3!
4SUBROUTINE STRATEMIT(pdtphys,pdt,xlat,xlon,t_seri,pplay,paprs,tr_seri,&
5                    m_emiss_vol_daily,latmin,latmax,lonmin,lonmax,altemiss,sigma_alt,id_spec,&
6                    stretchlong,ispeci,id_species_total)
7
8  USE dimphy, ONLY : klon,klev
9  USE strataer_local_var_mod
10  USE phys_cal_mod
11  USE phys_local_var_mod, ONLY: d_q_emiss
12  USE YOMCST, only : RD, RPI, RG
13  USE geometry_mod, ONLY : cell_area, boundslat
14  USE aerophys
15  USE infotrac_phy
16  USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
17  USE mod_grid_phy_lmdz, ONLY: nbp_lon
18 
19  IMPLICIT NONE
20
21! Input argument
22!---------------
23    REAL,INTENT(IN)    :: pdtphys, pdt               ! Pas d'integration pour la physique(seconde)
24    REAL,DIMENSION(klon),INTENT(IN)    :: xlat       ! latitudes pour chaque point
25    REAL,DIMENSION(klon),INTENT(IN)    :: xlon       ! longitudes pour chaque point
26    REAL,DIMENSION(klon,klev),INTENT(IN)   :: pplay
27    ! pression pour le mileu de chaque couche (en Pa)
28    REAL,DIMENSION(klon,klev+1),INTENT(IN) :: paprs
29    ! pression pour chaque inter-couche (en Pa)
30    REAL,DIMENSION(klon,klev),INTENT(IN)   :: t_seri ! Temperature
31    REAL,INTENT(IN)    :: latmin,latmax,lonmin,lonmax
32    REAL,INTENT(IN)    :: m_emiss_vol_daily          ! daily injection mass emission
33    REAL,INTENT(IN)    :: stretchlong                ! length to stretch emission (day)
34    REAL,DIMENSION(klon,klev) :: m_air_gridbox       ! mass ofair in every grid box[kg
35    REAL,INTENT(IN)     :: sigma_alt, altemiss
36    INTEGER,INTENT(IN)  :: id_spec, ispeci, id_species_total
37   
38! Output argument
39!----------------
40    REAL,DIMENSION(klon,klev,nbtr),INTENT(INOUT)  :: tr_seri ! Concentration Traceur [U/KgA]
41! Local variables
42!----------------
43    INTEGER              :: ieru, i, k, i_int, ilon, ilev
44    REAL                 :: dlat_loc
45    REAL,DIMENSION(klev+1)  :: altLMDz    ! altitude of layer interfaces in m
46    REAL                 :: zrho          ! Density of air [kg/m3]
47    REAL,DIMENSION(klev) :: zdm           ! mass of atm. model layer in kg 
48    REAL                 :: zdz           ! thickness of atm. model layer in m
49    REAL                 :: f_lay_sum     ! sum of layer emission fractions
50    REAL,DIMENSION(klev) :: f_lay_emiss   ! emiss fraction for every emiss for every vertical layer
51    REAL                 :: emission      ! emission
52    REAL                 :: theta_min, theta_max, theta ! for SAI computation between two latitudes
53   
54    IF (is_mpi_root) THEN
55       WRITE(*,*) 'IN STRATEMIT: date from phys_cal_mod=',year_cur,'-',&
56                      &  mth_cur,'-',day_cur,'-',hour,' flh2o=',flh2o
57    ENDIF
58
59!--calculate mass of air in every grid box !NL: DONE ABOVE DIRECTLY IN MAIN LOOP I,K
60!    DO ilon=1, klon
61!       DO ilev=1, klev
62!          m_air_gridbox(ilon,ilev)=(paprs(ilon,ilev)-paprs(ilon,ilev+1))/RG*cell_area(ilon)
63!       ENDDO
64!    ENDDO
65   
66    DO i=1,klon
67       
68       dlat_loc=180./RPI/2.*(boundslat(i,1)-boundslat(i,3))
69       dlon=360./2./FLOAT(nbp_lon)
70       ! dlat = half difference of boundary latitudes
71       IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT, dlat=',dlat_loc,xlat(i),xlon(i)
72       
73       theta_min = max(xlat(i)-dlat_loc,latmin)
74       theta_max = min(xlat(i)+dlat_loc,latmax)
75               
76       IF (  xlat(i).GE.latmin-dlat_loc .AND. &
77         &   xlat(i).LT.latmax+dlat_loc .AND. &
78         &   xlon(i).GE.lonmin-dlon .AND. &
79         &   xlon(i).LT.lonmax+dlon ) THEN
80                   !
81          WRITE(*,*) 'coordinates of volcanic injection point=',&
82            &    xlat(i),xlon(i),day_cur,mth_cur,year_cur
83          WRITE(*,*) 'DD m_emiss_vol_daily=', &
84            &    m_emiss_vol_daily
85         
86          !compute altLMDz
87          altLMDz(:)=0.0
88          DO k=1, klev
89             zrho=pplay(i,k)/t_seri(i,k)/RD       !air density in kg/m3
90             zdm(k)=(paprs(i,k)-paprs(i,k+1))/RG  !mass of layer in kg
91             zdz=zdm(k)/zrho                   !thickness of layer in m
92             altLMDz(k+1)=altLMDz(k)+zdz       !altitude of interface
93          ENDDO
94         
95          CALL STRATDISTRIB(altLMDz,altemiss,sigma_alt,f_lay_emiss)
96         
97          IF (flag_emit.EQ.3) then
98             theta=(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ &
99               & (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI))
100          ELSE
101             theta=1.
102          ENDIF
103
104          !vertically distributed emission
105          DO k=1, klev
106            !--calculate mass of air in every grid box
107            m_air_gridbox(i,k)=(paprs(i,k)-paprs(i,k+1))/RG*cell_area(i)
108            ! stretch emission over stretchlong period
109            emission=m_emiss_vol_daily/m_air_gridbox(i,k)*f_lay_emiss(k)/stretchlong/ &
110                 &     (86400.-pdt)*theta
111           
112             IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: emission avant/apres', &
113              &  'alt= ', altLMDz(k), &
114              &  'flh2o= ',flh2o, &
115              &  'id_speac= ',id_spec,id_species_total, &
116              &  'emission= ',emission, &
117              &  'pdtphys= ',pdtphys, &
118              &  'rapport m_emiss/m_air*f_lay= ', m_emiss_vol_daily/m_air_gridbox(i,k)*f_lay_emiss(k), &
119              &   'stretchlong= ', stretchlong, &
120              &   'theta= ', theta
121             
122             IF(emission < 1.E-34) emission = 0.0
123             
124             IF (flh2o.EQ.0) THEN
125                IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: tr_ser avant/apres',&
126                  &  'i= ',i,'k= ',k, 'flh2o= ',flh2o, &
127                  &  tr_seri(i,k,id_spec), &
128                  &  tr_seri(i,k,id_spec)+emission*pdtphys
129             
130                tr_seri(i,k,id_spec)=tr_seri(i,k,id_spec)+emission*pdtphys
131                IF (id_species_total.NE.0) THEN
132                   tr_seri(i,k,id_species_total)=tr_seri(i,k,id_species_total)+emission*pdtphys
133                ENDIF
134             ELSE IF(flh2o.EQ.1) THEN
135                d_q_emiss(i,k)=emission*pdtphys
136                IF(d_q_emiss(i,k) > 1.E34) THEN
137                   WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) > 1.E34 = ',d_q_emiss(i,k)
138                ELSE IF(d_q_emiss(i,k) < -1.E34) THEN
139                   WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) < -1.E34 = ',d_q_emiss(i,k)
140                ENDIF
141               
142                IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: ',&
143                     &  'i= ',i,'k= ',k, 'flh2o= ',flh2o, &
144                     &  'emission= ',emission, &
145                     &  'd_q_emiss(i,k)= ',d_q_emiss(i,k)
146               
147                IF(d_q_emiss(i,k) > 1.E34) THEN
148                   d_q_emiss(i,k) = 1.E34
149                   WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) > 1.E34 : ',d_q_emiss(i,k)
150                ELSE IF(d_q_emiss(i,k) < -1.E34) THEN
151                   d_q_emiss(i,k) = -1.E34
152                   WRITE(*,*) 'IN STRATEMIT: d_q_emiss(i,k) < -1.E34 : ',d_q_emiss(i,k)
153                ENDIF
154             ENDIF
155             budg_emi(i,ispeci)=budg_emi(i,ispeci)+emission*zdm(k)
156          ENDDO
157         
158       ENDIF ! emission grid cell
159
160    ENDDO ! klon loop
161
162END SUBROUTINE STRATEMIT
Note: See TracBrowser for help on using the repository browser.