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

Last change on this file since 1242 was 1226, checked in by aslmd, 11 years ago

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