source: LMDZ6/branches/Amaury_dev/libf/phylmd/StratAer/stratemit.F90 @ 5409

Last change on this file since 5409 was 5116, checked in by abarral, 5 months ago

rename modules properly lmdz_*
move ismin, ismax, minmax into new lmdz_libmath.f90
(lint) uppercase fortran keywords
File size: 6.7 KB
RevLine 
[5099]1
[4601]2! $Id: stratemit.F90  2022-07-04 mmarchand $
[5099]3
[4601]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
[5101]8  USE dimphy, ONLY: klon,klev
[4601]9  USE strataer_local_var_mod
10  USE phys_cal_mod
11  USE phys_local_var_mod, ONLY: d_q_emiss
[5101]12  USE lmdz_yomcst, ONLY: RD, RPI, RG
[5112]13  USE lmdz_geometry, ONLY: cell_area, boundslat
[4601]14  USE aerophys
15  USE infotrac_phy
[5110]16  USE lmdz_phys_mpi_data, ONLY:  is_mpi_root
17  USE lmdz_grid_phy, ONLY: nbp_lon
[4601]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,'-',&
[5087]56    mth_cur,'-',day_cur,'-',hour,' flh2o=',flh2o
[4601]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               
[5082]76       IF (  xlat(i)>=latmin-dlat_loc .AND. &
[5087]77     xlat(i)<latmax+dlat_loc .AND. &
78     xlon(i)>=lonmin-dlon .AND. &
79     xlon(i)<lonmax+dlon ) THEN
[5099]80
[4601]81          WRITE(*,*) 'coordinates of volcanic injection point=',&
[5087]82      xlat(i),xlon(i),day_cur,mth_cur,year_cur
[4601]83          WRITE(*,*) 'DD m_emiss_vol_daily=', &
[5087]84      m_emiss_vol_daily
[4601]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         
[5116]97          IF (flag_emit==3) THEN
[4601]98             theta=(sin(theta_max/180.*RPI)-sin(theta_min/180.*RPI))/ &
[5087]99   (sin(xlat_max_sai/180.*RPI)-sin(xlat_min_sai/180.*RPI))
[4601]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/ &
[5087]110       (86400.-pdt)*theta
[4601]111           
112             IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: emission avant/apres', &
[5087]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
[4601]121             
122             IF(emission < 1.E-34) emission = 0.0
123             
[5082]124             IF (flh2o==0) THEN
[4601]125                IF(flag_verbose_strataer) WRITE(*,*) 'IN STRATEMIT: tr_ser avant/apres',&
[5087]126    'i= ',i,'k= ',k, 'flh2o= ',flh2o, &
127    tr_seri(i,k,id_spec), &
128    tr_seri(i,k,id_spec)+emission*pdtphys
[4601]129             
130                tr_seri(i,k,id_spec)=tr_seri(i,k,id_spec)+emission*pdtphys
[5082]131                IF (id_species_total/=0) THEN
[4601]132                   tr_seri(i,k,id_species_total)=tr_seri(i,k,id_species_total)+emission*pdtphys
133                ENDIF
[5082]134             ELSE IF(flh2o==1) THEN
[4601]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: ',&
[5087]143    'i= ',i,'k= ',k, 'flh2o= ',flh2o, &
144    'emission= ',emission, &
145    'd_q_emiss(i,k)= ',d_q_emiss(i,k)
[4601]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.