source: trunk/LMDZ.MARS/libf/phymars/callsedim.F @ 1747

Last change on this file since 1747 was 1720, checked in by jaudouard, 7 years ago

Update on CO2 ice clouds scheme
File size: 16.7 KB
Line 
1      SUBROUTINE callsedim(ngrid,nlay, ptimestep,
2     &                pplev,zlev, zlay, pt, pdt, rdust, rice,
3     &                rsedcloud,rhocloud,
4     &                pq, pdqfi, pdqsed,pdqs_sed,nq,
5     &                tau,tauscaling)
6! to use  'getin'
7      USE ioipsl_getincom, only: getin
8      USE updaterad, only: updaterdust,updaterice_micro,updaterice_typ
9      USE tracer_mod, only: noms, igcm_dust_mass, igcm_dust_number,
10     &                      rho_dust, rho_q, radius, varian,
11     &                      igcm_ccn_mass, igcm_ccn_number,
12     &                      igcm_h2o_ice, nuice_sed, nuice_ref,
13     &                      igcm_ccnco2_mass, igcm_ccnco2_number,
14     &                      igcm_co2_ice
15      USE comcstfi_h
16      IMPLICIT NONE
17
18c=======================================================================
19c      Sedimentation of the  Martian aerosols
20c      depending on their density and radius
21c
22c      F.Forget 1999
23c
24c      Modified by J.-B. Madeleine 2010: Now includes the doubleq
25c        technique in order to have only one call to callsedim in
26c        physiq.F.
27c
28c      Modified by J. Audouard 09/16: Now includes the co2clouds case
29c        If the co2 microphysics is on, then co2 theice & ccn tracers
30c        are being sedimented in the microtimestep (co2cloud.F), not
31c        in this routine.
32c
33c=======================================================================
34
35c-----------------------------------------------------------------------
36c   declarations:
37c   -------------
38     
39#include "callkeys.h"
40
41c
42c   arguments:
43c   ----------
44
45      integer,intent(in) :: ngrid  ! number of horizontal grid points
46      integer,intent(in) :: nlay   ! number of atmospheric layers
47      real,intent(in) :: ptimestep ! physics time step (s)
48      real,intent(in) :: pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
49      real,intent(in) :: zlev(ngrid,nlay+1) ! altitude at layer boundaries
50      real,intent(in) :: zlay(ngrid,nlay)   ! altitude at the middle of the layers
51      real,intent(in) :: pt(ngrid,nlay) ! temperature at mid-layer (K)
52      real,intent(in) :: pdt(ngrid,nlay) ! tendency on temperature, from
53                                         ! previous processes (K/s)
54c    Aerosol radius provided by the water ice microphysical scheme:
55      real,intent(out) :: rdust(ngrid,nlay) ! Dust geometric mean radius (m)
56      real,intent(out) :: rice(ngrid,nlay)  ! H2O Ice geometric mean radius (m)
57c     Sedimentation radius of water ice
58      real,intent(in) :: rsedcloud(ngrid,nlay)
59c     Cloud density (kg.m-3)
60      real,intent(inout) :: rhocloud(ngrid,nlay)
61c    Traceurs :
62      real,intent(in) :: pq(ngrid,nlay,nq)  ! tracers (kg/kg)
63      real,intent(in) :: pdqfi(ngrid,nlay,nq)  ! tendency before sedimentation (kg/kg.s-1)
64      real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
65      real,intent(out) :: pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)
66      integer,intent(in) :: nq  ! number of tracers
67      real,intent(in) :: tau(ngrid,nlay) ! dust opacity
68      real,intent(in) :: tauscaling(ngrid)
69     
70c   local:
71c   ------
72
73      INTEGER l,ig, iq
74
75      real zqi(ngrid,nlay,nq) ! to locally store tracers
76      real zt(ngrid,nlay) ! to locally store temperature
77      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
78      real epaisseur (ngrid,nlay) ! Layer thickness (m)
79      real wq(ngrid,nlay+1) ! displaced tracer mass (kg.m-2)
80      real r0(ngrid,nlay) ! geometric mean radius used for
81                                !   sedimentation (m)
82      real r0dust(ngrid,nlay) ! geometric mean radius used for
83                                    !   dust (m)
84!      real r0ccn(ngrid,nlay)  ! geometric mean radius used for
85!                                    !   CCNs (m)
86      real,save :: beta ! correction for the shape of the ice particles (cf. newsedim)
87
88c     for ice radius computation
89      REAL Mo,No
90      REAl ccntyp
91
92
93
94c     Discrete size distributions (doubleq)
95c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
96c       1) Parameters used to represent the changes in fall
97c          velocity as a function of particle size;
98      integer ir
99      integer,parameter :: nr=12 !(nr=7) ! number of bins
100      real,save :: rd(nr)
101      real qr(ngrid,nlay,nr)
102      real,save :: rdi(nr+1)    ! extreme and intermediate radii
103      real Sq(ngrid,nlay)
104      real,parameter :: rdmin=1.e-8
105      real,parameter :: rdmax=30.e-6
106      real,parameter :: rdimin=1.e-8 ! 1.e-7
107      real,parameter :: rdimax=1.e-4
108
109c       2) Second size distribution for the log-normal integration
110c          (the mass mixing ratio is computed for each radius)
111
112      integer iint
113      integer,parameter :: ninter=4 ! number of points between each rdi radii
114      real,save :: rr(ninter,nr)
115      integer radpower
116      real sigma0
117
118c       3) Other local variables used in doubleq
119
120      INTEGER,SAVE :: idust_mass  ! index of tracer containing dust mass
121                                  !   mix. ratio
122      INTEGER,SAVE :: idust_number ! index of tracer containing dust number
123                                   !   mix. ratio
124      INTEGER,SAVE :: iccn_mass  ! index of tracer containing CCN mass
125                                 !   mix. ratio
126      INTEGER,SAVE :: iccn_number ! index of tracer containing CCN number
127                                  !   mix. ratio
128      INTEGER,SAVE :: iccnco2_number ! index of tracer containing CCN number
129      INTEGER,SAVE :: iccnco2_mass ! index of tracer containing CCN number
130      INTEGER,SAVE :: ico2_ice ! index of tracer containing CCN number
131
132
133      LOGICAL,SAVE :: firstcall=.true.
134
135c    ** un petit test de coherence
136c       --------------------------
137
138      IF (firstcall) THEN
139         
140c       Doubleq: initialization
141        IF (doubleq) THEN
142         do ir=1,nr
143             rd(ir)= rdmin*(rdmax/rdmin)**(float(ir-1)/float(nr-1))
144         end do
145         rdi(1)=rdimin
146         do ir=2,nr
147           rdi(ir)= sqrt(rd(ir-1)*rd(ir))
148         end do
149         rdi(nr+1)=rdimax
150
151         do ir=1,nr
152           do iint=1,ninter
153             rr(iint,ir)=
154     &        rdi(ir)*
155     &        (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
156c             write(*,*) rr(iint,ir)
157           end do
158         end do
159
160      ! identify tracers corresponding to mass mixing ratio and
161      ! number mixing ratio
162        idust_mass=0      ! dummy initialization
163        idust_number=0    ! dummy initialization
164
165        do iq=1,nq
166          if (noms(iq).eq."dust_mass") then
167            idust_mass=iq
168            write(*,*)"callsedim: idust_mass=",idust_mass
169          endif
170          if (noms(iq).eq."dust_number") then
171            idust_number=iq
172            write(*,*)"callsedim: idust_number=",idust_number
173          endif
174        enddo
175
176        ! check that we did find the tracers
177        if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
178          write(*,*) 'callsedim: error! could not identify'
179          write(*,*) ' tracers for dust mass and number mixing'
180          write(*,*) ' ratio and doubleq is activated!'
181          stop
182        endif
183        ENDIF !of if (doubleq)
184
185        IF (microphys) THEN
186          iccn_mass=0
187          iccn_number=0
188          do iq=1,nq
189            if (noms(iq).eq."ccn_mass") then
190              iccn_mass=iq
191              write(*,*)"callsedim: iccn_mass=",iccn_mass
192            endif
193            if (noms(iq).eq."ccn_number") then
194              iccn_number=iq
195              write(*,*)"callsedim: iccn_number=",iccn_number
196            endif
197          enddo
198          ! check that we did find the tracers
199          if ((iccn_mass.eq.0).or.(iccn_number.eq.0)) then
200            write(*,*) 'callsedim: error! could not identify'
201            write(*,*) ' tracers for ccn mass and number mixing'
202            write(*,*) ' ratio and microphys is activated!'
203            stop
204          endif
205        ENDIF !of if (microphys)
206
207        IF (co2clouds) THEN
208          iccnco2_mass=0
209          iccnco2_number=0
210          ico2_ice=0
211          do iq=1,nq
212            if (noms(iq).eq."ccnco2_mass") then
213              iccnco2_mass=iq
214              write(*,*)"callsedim: iccnco2_mass=",iccnco2_mass
215            endif
216            if (noms(iq).eq."co2_ice") then
217              ico2_ice=iq
218              write(*,*)"callsedim: ico2_ice=",ico2_ice
219            endif
220            if (noms(iq).eq."ccnco2_number") then
221              iccnco2_number=iq
222              write(*,*)"callsedim: iccnco2_number=",iccnco2_number
223            endif
224          enddo
225          ! check that we did find the tracers
226          if ((iccnco2_mass.eq.0).or.(iccnco2_number.eq.0)) then
227            write(*,*) 'callsedim: error! could not identify'
228            write(*,*) ' tracers for ccn co2 mass and number mixing'
229            write(*,*) ' ratio and co2clouds are activated!'
230            stop
231          endif
232       ENDIF                    !of if (co2clouds)
233
234
235        IF (water) THEN
236         write(*,*) "correction for the shape of the ice particles ?"
237         beta=0.75 ! default value
238         call getin("ice_shape",beta)
239         write(*,*) " ice_shape = ",beta
240
241          write(*,*) "water_param nueff Sedimentation:", nuice_sed
242          IF (activice) THEN
243            write(*,*) "water_param nueff Radiative:", nuice_ref
244          ENDIF
245        ENDIF
246     
247        firstcall=.false.
248      ENDIF ! of IF (firstcall)
249
250c-----------------------------------------------------------------------
251c    1. Initialization
252c    -----------------
253
254!      zqi(1:ngrid,1:nlay,1:nqmx) = 0.
255c     Update the mass mixing ratio and temperature with the tendencies coming
256c       from other parameterizations:
257c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
258      zqi(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)
259     &                         +pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
260      zt(1:ngrid,1:nlay)=pt(1:ngrid,1:nlay)
261     &                         +pdt(1:ngrid,1:nlay)*ptimestep
262
263
264c    Computing the different layer properties
265c    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
266c    Mass (kg.m-2), thickness(m), crossing time (s)  etc.
267
268      do  l=1,nlay
269        do ig=1, ngrid
270          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
271          epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
272        end do
273      end do
274
275c =================================================================
276c     Compute the geometric mean radius used for sedimentation
277
278      if (doubleq) then
279        do l=1,nlay
280          do ig=1, ngrid
281     
282         call updaterdust(zqi(ig,l,igcm_dust_mass),
283     &                    zqi(ig,l,igcm_dust_number),r0dust(ig,l),
284     &                    tauscaling(ig))
285         
286          end do
287        end do
288      endif
289
290
291c =================================================================
292      do iq=1,nq
293        if(radius(iq).gt.1.e-9 .and.(iq.ne.ico2_ice) .and.
294     &        (iq .ne. iccnco2_mass) .and. (iq .ne.
295     &        iccnco2_number)) then   ! no sedim for gaz or CO2 clouds  (done in microtimestep)
296
297c -----------------------------------------------------------------
298c         DOUBLEQ CASE
299c -----------------------------------------------------------------
300          if ((doubleq.and.
301     &        ((iq.eq.idust_mass).or.
302     &         (iq.eq.idust_number)))) then
303     
304c           Computing size distribution:
305c           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
306
307c            if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
308              do  l=1,nlay
309                do ig=1, ngrid
310                  r0(ig,l)=r0dust(ig,l)
311                end do
312              end do
313              sigma0 = varian
314
315c        Computing mass mixing ratio for each particle size
316c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
317          IF ((iq.EQ.idust_mass).or.(iq.EQ.iccn_mass)) then
318            radpower = 2
319          ELSE  ! number
320            radpower = -1
321          ENDIF
322          Sq(1:ngrid,1:nlay) = 0.
323          do ir=1,nr
324            do l=1,nlay
325              do ig=1,ngrid
326c                ****************
327c                Size distribution integration
328c                (Trapezoid Integration Method)
329                 qr(ig,l,ir)=0.5*(rr(2,ir)-rr(1,ir))*
330     &             (rr(1,ir)**radpower)*
331     &             exp(-(log(rr(1,ir)/r0(ig,l)))**2/(2*sigma0**2))
332                 do iint=2,ninter-1
333                   qr(ig,l,ir)=qr(ig,l,ir) +
334     &             0.5*(rr(iint+1,ir)-rr(iint-1,ir))*
335     &             (rr(iint,ir)**radpower)*
336     &             exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
337     &             (2*sigma0**2))
338                 end do
339                 qr(ig,l,ir)=qr(ig,l,ir) +
340     &             0.5*(rr(ninter,ir)-rr(ninter-1,ir))*
341     &             (rr(ninter,ir)**radpower)*
342     &             exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
343     &             (2*sigma0**2))
344
345c                **************** old method (not recommended!)
346c                qr(ig,l,ir)=(rd(ir)**(5-3*iq))*
347c    &           exp( -(log(rd(ir)/r0(ig,l)))**2 / (2*sigma0**2) )
348c                ******************************
349
350                 Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
351              enddo
352            enddo
353          enddo
354
355          do ir=1,nr
356            do l=1,nlay
357              do ig=1,ngrid
358                 qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
359              enddo
360            enddo
361          enddo
362
363c         Computing sedimentation for each tracer
364c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
365
366          zqi(1:ngrid,1:nlay,iq) = 0.
367          pdqs_sed(1:ngrid,iq) = 0.
368
369          do ir=1,nr
370         
371               call newsedim(ngrid,nlay,1,1,ptimestep,
372     &         pplev,masse,epaisseur,zt,rd(ir),rho_dust,qr(1,1,ir),
373     &         wq,0.5)
374
375c            Tendencies
376c            ~~~~~~~~~~
377             do ig=1,ngrid
378               pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
379     &                                + wq(ig,1)/ptimestep
380             end do
381             DO l = 1, nlay
382               DO ig=1,ngrid
383                 zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
384               ENDDO
385             ENDDO
386          enddo ! of do ir=1,nr
387c -----------------------------------------------------------------
388c         WATER CYCLE CASE
389c -----------------------------------------------------------------
390           else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
391     &       .or. (iq .eq. igcm_h2o_ice)) then
392            if (microphys) then
393              ! water ice sedimentation
394              call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
395     &        ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rhocloud,
396     &        zqi(1,1,iq),wq,beta)
397            else
398              ! water ice sedimentation
399              call newsedim(ngrid,nlay,ngrid*nlay,1,
400     &        ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rho_q(iq),
401     &        zqi(1,1,iq),wq,beta)
402            endif ! of if (microphys)
403c           Tendencies
404c           ~~~~~~~~~~
405            do ig=1,ngrid
406              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
407            end do
408c -----------------------------------------------------------------
409c         GENERAL CASE
410c -----------------------------------------------------------------
411          else
412            call newsedim(ngrid,nlay,1,1,ptimestep,
413     &      pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
414     &      zqi(1,1,iq),wq,1.0)
415c           Tendencies
416c           ~~~~~~~~~~
417            do ig=1,ngrid
418              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
419            end do
420          endif ! of if doubleq and if water
421c -----------------------------------------------------------------
422
423c         Compute the final tendency:
424c         ---------------------------
425          DO l = 1, nlay
426            DO ig=1,ngrid
427              pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
428     $        (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
429            ENDDO
430          ENDDO
431
432        endif ! of if(radius(iq).gt.1.e-9)
433c =================================================================
434      enddo ! of do iq=1,nq
435
436c     Update the dust particle size "rdust"
437c     -------------------------------------
438      if (doubleq) then
439       DO l = 1, nlay
440        DO ig=1,ngrid
441       
442     
443         call updaterdust(zqi(ig,l,igcm_dust_mass),
444     &                    zqi(ig,l,igcm_dust_number),rdust(ig,l),
445     &                    tauscaling(ig))     
446
447         
448        ENDDO
449       ENDDO
450      endif ! of if (doubleq)
451     
452c     Update the ice particle size "rice"
453c     -------------------------------------
454      if (water) then
455       IF(microphys) THEN
456       
457       
458        DO l = 1, nlay
459          DO ig=1,ngrid
460
461         call updaterice_micro(zqi(ig,l,igcm_h2o_ice),
462     &    zqi(ig,l,igcm_ccn_mass),zqi(ig,l,igcm_ccn_number),
463     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
464           
465          ENDDO
466        ENDDO
467       
468       ELSE
469       
470        DO l = 1, nlay
471          DO ig=1,ngrid
472         
473            call updaterice_typ(zqi(ig,l,igcm_h2o_ice),
474     &                      tau(ig,1),zlay(ig,l),rice(ig,l))
475
476          ENDDO
477        ENDDO
478       ENDIF ! of IF(microphys)
479      endif ! of if (water)
480
481      END
482
Note: See TracBrowser for help on using the repository browser.