source: trunk/LMDZ.VENUS/libf/phyvenus/phytrac_chimie.F @ 2487

Last change on this file since 2487 was 2464, checked in by slebonnois, 4 years ago

SL: major update related to extension to 250 km. New EUV heating as used in Martian GCM, debug of a few routines, adding He, clean up, updating mmean regularly, modification of the order of processes, add the possibility to use Hedin composition in rad transfer

  • Property svn:executable set to *
File size: 11.3 KB
Line 
1      subroutine phytrac_chimie (
2     $                    debutphy,
3     $                    gmtime,
4     $                    nqmax,
5     $                    nlon,
6     $                    lat,
7     $                    lon,
8     $                    nlev,
9     $                    pdtphys,
10     $                    temp,
11     $                    pplay,
12     $                    trac,
13     $                    d_tr_chem,
14     $                    iter)
15
16      use chemparam_mod
17      use conc, only: mmean,rnew
18
19      implicit none
20     
21#include "clesphys.h"
22#include "YOMCST.h"
23
24!===================================================================
25!     input
26!===================================================================
27
28      integer :: nlon, nlev                     ! number of gridpoints and levels
29      integer :: nqmax                          ! number of tracers
30
31      real :: gmtime                            ! day fraction
32      real :: pdtphys                           ! phytrac_chimie timestep (s)
33      real, dimension(nlon,nlev) :: temp        ! temperature (k)
34      real, dimension(nlon,nlev) :: pplay       ! pressure (pa)
35      real, dimension(nlon,nlev,nqmax) :: trac  ! tracer mass mixing ratio
36
37      logical :: debutphy                       ! first call flag
38
39!===================================================================
40!     output
41!===================================================================
42
43      real, dimension(nlon,nlev,nqmax) :: d_tr_chem  ! chemical tendency for each tracer
44      integer, dimension(nlon,nlev) :: iter          ! chemical iterations
45
46!===================================================================
47!     local
48!===================================================================
49
50      real :: sza_local   ! solar zenith angle (deg)
51      real :: lon_sun
52
53      integer :: i, iq
54      integer :: ilon, ilev
55
56      real  lat(nlon), lat_local(nlon)
57      real  lon(nlon), lon_local(nlon)
58
59      real, dimension(nlon,nlev) :: mrtwv, mrtsa ! total water and total sulfuric acid
60      real, dimension(nlon,nlev) :: mrwv, mrsa   ! gas-phase water and gas-phase sulfuric acid
61      real, dimension(nlon,nlev) :: trac_sum
62      real, dimension(nlon,nlev,nqmax) :: ztrac  ! local tracer mixing ratio
63     
64!===================================================================
65!     first call : initialisations
66!===================================================================
67
68      if (debutphy) then
69     
70!-------------------------------------------------------------------
71!        case of tracers re-initialisation with chemistry
72!-------------------------------------------------------------------
73         if (reinit_trac .and. ok_chem) then
74
75!!! in this reinitialisation, trac is VOLUME mixing ratio
76            print*, "Tracers are re-initialised"
77            trac(:,:,:) = 1.0e-30
78
79            if ((i_ocs /= 0) .and. (i_co /= 0) .and. (i_hcl /= 0)
80     $           .and. (i_so2 /= 0) .and. (i_h2o /= 0) .and. (i_n2/ = 0)
81     $           .and. (i_co2 /= 0)) then
82
83               trac(:,1:22,i_ocs) = 3.e-6
84               trac(:,1:22,i_co)  = 25.e-6
85               trac(:,:,i_hcl)    = 0.4e-6
86               trac(:,1:22,i_so2) = 7.e-6
87               trac(:,1:22,i_h2o) = 30.e-6
88               trac(:,:,i_n2)     = 0.35e-1
89   
90!          ensure that sum of mixing ratios = 1
91
92               trac_sum(:,:) = 0.
93
94               do iq = 1,nqmax - nmicro
95                  if (iq /= i_co2) then
96                     trac_sum(:,:) = trac_sum(:,:) + trac(:,:,iq)
97                  end if
98               end do
99
100!          initialise co2
101
102               trac(:,:,i_co2) = 1. - trac_sum(:,:)
103
104            else
105               write(*,*) "at least one tracer is missing : stop"
106               stop
107            end if
108       
109
110! update mmean
111            mmean(:,:) = 0.
112            do iq = 1,nqmax - nmicro
113               mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
114               rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
115            enddo
116
117!           convert volume to mass mixing ratio
118            do iq = 1,nqmax - nmicro
119               trac(:,:,iq) = trac(:,:,iq)*m_tr(iq)/mmean(:,:)
120            end do
121   
122         end if  ! reinit_trac
123
124!-------------------------------------------------------------------
125!        case of detailed microphysics without chemistry
126!-------------------------------------------------------------------
127         if (.not. ok_chem .and. ok_cloud .and. cl_scheme == 2) then
128
129!           convert mass to volume mixing ratio
130
131            do iq = 1,nqmax - nmicro
132               ztrac(:,:,iq) = trac(:,:,iq)*mmean(:,:)/m_tr(iq)
133            end do
134         
135!           initialise microphysics
136 
137            call vapors4muphy_ini(nlon,nlev,ztrac)
138
139!           convert volume to mass mixing ratio
140
141            do iq = 1,nqmax - nmicro
142               trac(:,:,iq) = ztrac(:,:,iq)*m_tr(iq)/mmean(:,:)
143            end do
144   
145         end if
146
147      end if  ! debutphy
148
149!===================================================================
150!     convert mass to volume mixing ratio : gas phase
151!===================================================================
152
153      do iq = 1,nqmax - nmicro
154         ztrac(:,:,iq) = max(trac(:,:,iq)*mmean(:,:)/m_tr(iq), 1.e-30)
155      end do
156
157!===================================================================
158!     microphysics: simplified scheme (phd aurelien stolzenbach)
159!===================================================================
160
161      if (ok_cloud .and. cl_scheme == 1) then
162
163!     convert mass to volume mixing ratio : liquid phase
164
165         ztrac(:,:,i_h2so4liq) = max(trac(:,:,i_h2so4liq)
166     $                             *mmean(:,:)/m_tr(i_h2so4liq), 1.e-30)
167         ztrac(:,:,i_h2oliq) = max(trac(:,:,i_h2oliq)
168     $                             *mmean(:,:)/m_tr(i_h2oliq), 1.e-30)
169             
170!     total water and sulfuric acid (gas + liquid)
171
172         mrtwv(:,:) = ztrac(:,:,i_h2o) + ztrac(:,:,i_h2oliq)
173         mrtsa(:,:) = ztrac(:,:,i_h2so4) + ztrac(:,:,i_h2so4liq)
174
175!     all water and sulfuric acid is put in the gas-phase
176
177         mrwv(:,:) = mrtwv(:,:)
178         mrsa(:,:) = mrtsa(:,:)
179
180!     call microphysics
181
182         call new_cloud_venus(nlev, nlon, temp, pplay,
183     $                        mrtwv, mrtsa, mrwv, mrsa)
184
185!     update water vapour and sulfuric acid
186
187         ztrac(:,:,i_h2o) = mrwv(:,:)
188         ztrac(:,:,i_h2oliq) = mrtwv(:,:) - ztrac(:,:,i_h2o)
189     
190         ztrac(:,:,i_h2so4) = mrsa(:,:)
191         ztrac(:,:,i_h2so4liq) = mrtsa(:,:) - ztrac(:,:,i_h2so4)
192
193      end if  ! simplified scheme
194
195!===================================================================
196!     microphysics: detailed scheme (phd sabrina guilbon)
197!     !!! to be confirmed whether mad_muphy expects mmr or vmr for h2o and h2so4
198!===================================================================
199
200      if (ok_cloud .and. cl_scheme == 2) then
201
202         do iq = nqmax-nmicro+1,nqmax
203            ztrac(:,:,iq) = trac(:,:,iq)
204         end do
205
206         do ilon = 1,nlon       
207            do ilev = 1, nlev
208               if (temp(ilon,ilev) < 500.) then
209                  call mad_muphy(pdtphys,                               ! timestep
210     $                           temp(ilon,ilev),pplay(ilon,ilev),      ! temperature and pressure
211     $                           ztrac(ilon,ilev,i_h2o),
212     $                           ztrac(ilon,ilev,i_h2so4),     
213     $                           ztrac(ilon,ilev,i_m0_aer),
214     $                           ztrac(ilon,ilev,i_m3_aer),   
215     $                           ztrac(ilon,ilev,i_m0_mode1drop),
216     $                           ztrac(ilon,ilev,i_m0_mode1ccn),
217     $                           ztrac(ilon,ilev,i_m3_mode1sa),
218     $                           ztrac(ilon,ilev,i_m3_mode1w),   
219     $                           ztrac(ilon,ilev,i_m3_mode1ccn),   
220     $                           ztrac(ilon,ilev,i_m0_mode2drop),
221     $                           ztrac(ilon,ilev,i_m0_mode2ccn),
222     $                           ztrac(ilon,ilev,i_m3_mode2sa),
223     $                           ztrac(ilon,ilev,i_m3_mode2w),
224     $                           ztrac(ilon,ilev,i_m3_mode2ccn))
225               else
226                  ztrac(ilon,ilev,i_m0_aer)       = 0.
227                  ztrac(ilon,ilev,i_m3_aer)       = 0.
228                  ztrac(ilon,ilev,i_m0_mode1drop) = 0.
229                  ztrac(ilon,ilev,i_m0_mode1ccn)  = 0.
230                  ztrac(ilon,ilev,i_m3_mode1sa)   = 0.
231                  ztrac(ilon,ilev,i_m3_mode1w)    = 0.
232                  ztrac(ilon,ilev,i_m3_mode1ccn)  = 0.
233                  ztrac(ilon,ilev,i_m0_mode2drop) = 0.
234                  ztrac(ilon,ilev,i_m0_mode2ccn)  = 0.
235                  ztrac(ilon,ilev,i_m3_mode2sa)   = 0.
236                  ztrac(ilon,ilev,i_m3_mode2w)    = 0.
237                  ztrac(ilon,ilev,i_m3_mode2ccn)  = 0.
238               end if
239            end do
240         end do
241
242      end if  ! detailed scheme
243           
244!===================================================================
245!     photochemistry
246!===================================================================
247
248      if (ok_chem) then
249
250         lon_sun = (0.5 - gmtime)*2.*rpi
251         lon_local(:) = lon(:)*rpi/180.
252         lat_local(:) = lat(:)*rpi/180.
253
254         do ilon = 1,nlon
255
256!           solar zenith angle
257
258            sza_local = acos(cos(lat_local(ilon))*cos(lon_local(ilon))
259     $                 *cos(lon_sun) + cos(lat_local(ilon))
260     $                 *sin(lon_local(ilon))*sin(lon_sun))*180./rpi
261     
262            call photochemistry_venus(nlev, nlon, pdtphys,
263     $                                pplay(ilon,:)/100.,
264     $                                temp(ilon,:),
265     $                                ztrac(ilon,:,:),
266     $                                mmean(ilon,:),
267     $                                sza_local,
268     $                                lon(ilon), lat(ilon),
269     $                                nqmax, iter(ilon,:))
270
271         end do
272
273      end if  ! ok_chem
274
275!===================================================================
276!     compute tendencies
277!===================================================================
278
279! update mmean
280            mmean(:,:) = 0.
281            do iq = 1,nqmax - nmicro
282               mmean(:,:) = mmean(:,:)+trac(:,:,iq)*m_tr(iq)
283               rnew(:,:) = 8.314/mmean(:,:)*1.e3     ! J/kg K
284            enddo
285
286!     gas phase
287
288      do iq = 1,nqmax - nmicro
289         ztrac(:,:,iq) = max(ztrac(:,:,iq)*m_tr(iq)/mmean(:,:),
290     $                       1.e-30)
291         d_tr_chem(:,:,iq) = (ztrac(:,:,iq) - trac(:,:,iq))/pdtphys
292      end do
293
294!     liquid phase or moments
295
296      if (ok_cloud .and. cl_scheme == 1) then
297         ztrac(:,:,i_h2so4liq) = max(ztrac(:,:,i_h2so4liq)
298     $                               *m_tr(i_h2so4liq)/mmean(:,:),
299     $                               1.e-30)
300         ztrac(:,:,i_h2oliq)   = max(ztrac(:,:,i_h2oliq)
301     $                               *m_tr(i_h2oliq)/mmean(:,:),
302     $                               1.e-30)
303         d_tr_chem(:,:,i_h2so4liq) = (ztrac(:,:,i_h2so4liq)
304     $                              - trac(:,:,i_h2so4liq))/pdtphys
305         d_tr_chem(:,:,i_h2oliq) = (ztrac(:,:,i_h2oliq)
306     $                            - trac(:,:,i_h2oliq))/pdtphys
307      end if
308
309
310      if (ok_cloud .and. cl_scheme == 2) then
311         do iq = nqmax-nmicro+1,nqmax
312            d_tr_chem(:,:,iq) = (ztrac(:,:,iq) - trac(:,:,iq))/pdtphys
313         end do
314      end if
315
316      end subroutine phytrac_chimie
Note: See TracBrowser for help on using the repository browser.