source: trunk/LMDZ.MARS/libf/phymars/callsedim_mod.F @ 1962

Last change on this file since 1962 was 1962, checked in by mvals, 6 years ago

Mars GCM:

  • callsedim.F changed to module callsedim_mod.F

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