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

Last change on this file since 1818 was 1779, checked in by aslmd, 7 years ago

LMDZ.MARS (purely comments) marked the absolute firstcalls not supposed to change with runtime (e.g. not domain-related). this is most of them. those firstcall can stay local and do not need to be linked with the caller's general firstcall.
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      ! AS: firstcall OK absolute
139      IF (firstcall) THEN
140         
141c       Doubleq: initialization
142        IF (doubleq) THEN
143         do ir=1,nr
144             rd(ir)= rdmin*(rdmax/rdmin)**(float(ir-1)/float(nr-1))
145         end do
146         rdi(1)=rdimin
147         do ir=2,nr
148           rdi(ir)= sqrt(rd(ir-1)*rd(ir))
149         end do
150         rdi(nr+1)=rdimax
151
152         do ir=1,nr
153           do iint=1,ninter
154             rr(iint,ir)=
155     &        rdi(ir)*
156     &        (rdi(ir+1)/rdi(ir))**(float(iint-1)/float(ninter-1))
157c             write(*,*) rr(iint,ir)
158           end do
159         end do
160
161      ! identify tracers corresponding to mass mixing ratio and
162      ! number mixing ratio
163        idust_mass=0      ! dummy initialization
164        idust_number=0    ! dummy initialization
165
166        do iq=1,nq
167          if (noms(iq).eq."dust_mass") then
168            idust_mass=iq
169            write(*,*)"callsedim: idust_mass=",idust_mass
170          endif
171          if (noms(iq).eq."dust_number") then
172            idust_number=iq
173            write(*,*)"callsedim: idust_number=",idust_number
174          endif
175        enddo
176
177        ! check that we did find the tracers
178        if ((idust_mass.eq.0).or.(idust_number.eq.0)) then
179          write(*,*) 'callsedim: error! could not identify'
180          write(*,*) ' tracers for dust mass and number mixing'
181          write(*,*) ' ratio and doubleq is activated!'
182          stop
183        endif
184        ENDIF !of if (doubleq)
185
186        IF (microphys) THEN
187          iccn_mass=0
188          iccn_number=0
189          do iq=1,nq
190            if (noms(iq).eq."ccn_mass") then
191              iccn_mass=iq
192              write(*,*)"callsedim: iccn_mass=",iccn_mass
193            endif
194            if (noms(iq).eq."ccn_number") then
195              iccn_number=iq
196              write(*,*)"callsedim: iccn_number=",iccn_number
197            endif
198          enddo
199          ! check that we did find the tracers
200          if ((iccn_mass.eq.0).or.(iccn_number.eq.0)) then
201            write(*,*) 'callsedim: error! could not identify'
202            write(*,*) ' tracers for ccn mass and number mixing'
203            write(*,*) ' ratio and microphys is activated!'
204            stop
205          endif
206        ENDIF !of if (microphys)
207
208        IF (co2clouds) THEN
209          iccnco2_mass=0
210          iccnco2_number=0
211          ico2_ice=0
212          do iq=1,nq
213            if (noms(iq).eq."ccnco2_mass") then
214              iccnco2_mass=iq
215              write(*,*)"callsedim: iccnco2_mass=",iccnco2_mass
216            endif
217            if (noms(iq).eq."co2_ice") then
218              ico2_ice=iq
219              write(*,*)"callsedim: ico2_ice=",ico2_ice
220            endif
221            if (noms(iq).eq."ccnco2_number") then
222              iccnco2_number=iq
223              write(*,*)"callsedim: iccnco2_number=",iccnco2_number
224            endif
225          enddo
226          ! check that we did find the tracers
227          if ((iccnco2_mass.eq.0).or.(iccnco2_number.eq.0)) then
228            write(*,*) 'callsedim: error! could not identify'
229            write(*,*) ' tracers for ccn co2 mass and number mixing'
230            write(*,*) ' ratio and co2clouds are activated!'
231            stop
232          endif
233       ENDIF                    !of if (co2clouds)
234
235
236        IF (water) THEN
237         write(*,*) "correction for the shape of the ice particles ?"
238         beta=0.75 ! default value
239         call getin("ice_shape",beta)
240         write(*,*) " ice_shape = ",beta
241
242          write(*,*) "water_param nueff Sedimentation:", nuice_sed
243          IF (activice) THEN
244            write(*,*) "water_param nueff Radiative:", nuice_ref
245          ENDIF
246        ENDIF
247     
248        firstcall=.false.
249      ENDIF ! of IF (firstcall)
250
251c-----------------------------------------------------------------------
252c    1. Initialization
253c    -----------------
254
255!      zqi(1:ngrid,1:nlay,1:nqmx) = 0.
256c     Update the mass mixing ratio and temperature with the tendencies coming
257c       from other parameterizations:
258c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
259      zqi(1:ngrid,1:nlay,1:nq)=pq(1:ngrid,1:nlay,1:nq)
260     &                         +pdqfi(1:ngrid,1:nlay,1:nq)*ptimestep
261      zt(1:ngrid,1:nlay)=pt(1:ngrid,1:nlay)
262     &                         +pdt(1:ngrid,1:nlay)*ptimestep
263
264
265c    Computing the different layer properties
266c    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
267c    Mass (kg.m-2), thickness(m), crossing time (s)  etc.
268
269      do  l=1,nlay
270        do ig=1, ngrid
271          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
272          epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
273        end do
274      end do
275
276c =================================================================
277c     Compute the geometric mean radius used for sedimentation
278
279      if (doubleq) then
280        do l=1,nlay
281          do ig=1, ngrid
282     
283         call updaterdust(zqi(ig,l,igcm_dust_mass),
284     &                    zqi(ig,l,igcm_dust_number),r0dust(ig,l),
285     &                    tauscaling(ig))
286         
287          end do
288        end do
289      endif
290
291
292c =================================================================
293      do iq=1,nq
294        if(radius(iq).gt.1.e-9 .and.(iq.ne.ico2_ice) .and.
295     &        (iq .ne. iccnco2_mass) .and. (iq .ne.
296     &        iccnco2_number)) then   ! no sedim for gaz or CO2 clouds  (done in microtimestep)
297
298c -----------------------------------------------------------------
299c         DOUBLEQ CASE
300c -----------------------------------------------------------------
301          if ((doubleq.and.
302     &        ((iq.eq.idust_mass).or.
303     &         (iq.eq.idust_number)))) then
304     
305c           Computing size distribution:
306c           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
307
308c            if ((iq.eq.idust_mass).or.(iq.eq.idust_number)) then
309              do  l=1,nlay
310                do ig=1, ngrid
311                  r0(ig,l)=r0dust(ig,l)
312                end do
313              end do
314              sigma0 = varian
315
316c        Computing mass mixing ratio for each particle size
317c        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
318          IF ((iq.EQ.idust_mass).or.(iq.EQ.iccn_mass)) then
319            radpower = 2
320          ELSE  ! number
321            radpower = -1
322          ENDIF
323          Sq(1:ngrid,1:nlay) = 0.
324          do ir=1,nr
325            do l=1,nlay
326              do ig=1,ngrid
327c                ****************
328c                Size distribution integration
329c                (Trapezoid Integration Method)
330                 qr(ig,l,ir)=0.5*(rr(2,ir)-rr(1,ir))*
331     &             (rr(1,ir)**radpower)*
332     &             exp(-(log(rr(1,ir)/r0(ig,l)))**2/(2*sigma0**2))
333                 do iint=2,ninter-1
334                   qr(ig,l,ir)=qr(ig,l,ir) +
335     &             0.5*(rr(iint+1,ir)-rr(iint-1,ir))*
336     &             (rr(iint,ir)**radpower)*
337     &             exp(-(log(rr(iint,ir)/r0(ig,l)))**2/
338     &             (2*sigma0**2))
339                 end do
340                 qr(ig,l,ir)=qr(ig,l,ir) +
341     &             0.5*(rr(ninter,ir)-rr(ninter-1,ir))*
342     &             (rr(ninter,ir)**radpower)*
343     &             exp(-(log(rr(ninter,ir)/r0(ig,l)))**2/
344     &             (2*sigma0**2))
345
346c                **************** old method (not recommended!)
347c                qr(ig,l,ir)=(rd(ir)**(5-3*iq))*
348c    &           exp( -(log(rd(ir)/r0(ig,l)))**2 / (2*sigma0**2) )
349c                ******************************
350
351                 Sq(ig,l)=Sq(ig,l)+qr(ig,l,ir)
352              enddo
353            enddo
354          enddo
355
356          do ir=1,nr
357            do l=1,nlay
358              do ig=1,ngrid
359                 qr(ig,l,ir) = zqi(ig,l,iq)*qr(ig,l,ir)/Sq(ig,l)
360              enddo
361            enddo
362          enddo
363
364c         Computing sedimentation for each tracer
365c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
366
367          zqi(1:ngrid,1:nlay,iq) = 0.
368          pdqs_sed(1:ngrid,iq) = 0.
369
370          do ir=1,nr
371         
372               call newsedim(ngrid,nlay,1,1,ptimestep,
373     &         pplev,masse,epaisseur,zt,rd(ir),rho_dust,qr(1,1,ir),
374     &         wq,0.5)
375
376c            Tendencies
377c            ~~~~~~~~~~
378             do ig=1,ngrid
379               pdqs_sed(ig,iq) = pdqs_sed(ig,iq)
380     &                                + wq(ig,1)/ptimestep
381             end do
382             DO l = 1, nlay
383               DO ig=1,ngrid
384                 zqi(ig,l,iq)=zqi(ig,l,iq)+qr(ig,l,ir)
385               ENDDO
386             ENDDO
387          enddo ! of do ir=1,nr
388c -----------------------------------------------------------------
389c         WATER CYCLE CASE
390c -----------------------------------------------------------------
391           else if ((iq .eq. iccn_mass) .or. (iq .eq. iccn_number)
392     &       .or. (iq .eq. igcm_h2o_ice)) then
393            if (microphys) then
394              ! water ice sedimentation
395              call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
396     &        ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rhocloud,
397     &        zqi(1,1,iq),wq,beta)
398            else
399              ! water ice sedimentation
400              call newsedim(ngrid,nlay,ngrid*nlay,1,
401     &        ptimestep,pplev,masse,epaisseur,zt,rsedcloud,rho_q(iq),
402     &        zqi(1,1,iq),wq,beta)
403            endif ! of if (microphys)
404c           Tendencies
405c           ~~~~~~~~~~
406            do ig=1,ngrid
407              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
408            end do
409c -----------------------------------------------------------------
410c         GENERAL CASE
411c -----------------------------------------------------------------
412          else
413            call newsedim(ngrid,nlay,1,1,ptimestep,
414     &      pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
415     &      zqi(1,1,iq),wq,1.0)
416c           Tendencies
417c           ~~~~~~~~~~
418            do ig=1,ngrid
419              pdqs_sed(ig,iq)=wq(ig,1)/ptimestep
420            end do
421          endif ! of if doubleq and if water
422c -----------------------------------------------------------------
423
424c         Compute the final tendency:
425c         ---------------------------
426          DO l = 1, nlay
427            DO ig=1,ngrid
428              pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
429     $        (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
430            ENDDO
431          ENDDO
432
433        endif ! of if(radius(iq).gt.1.e-9)
434c =================================================================
435      enddo ! of do iq=1,nq
436
437c     Update the dust particle size "rdust"
438c     -------------------------------------
439      if (doubleq) then
440       DO l = 1, nlay
441        DO ig=1,ngrid
442       
443     
444         call updaterdust(zqi(ig,l,igcm_dust_mass),
445     &                    zqi(ig,l,igcm_dust_number),rdust(ig,l),
446     &                    tauscaling(ig))     
447
448         
449        ENDDO
450       ENDDO
451      endif ! of if (doubleq)
452     
453c     Update the ice particle size "rice"
454c     -------------------------------------
455      if (water) then
456       IF(microphys) THEN
457       
458       
459        DO l = 1, nlay
460          DO ig=1,ngrid
461
462         call updaterice_micro(zqi(ig,l,igcm_h2o_ice),
463     &    zqi(ig,l,igcm_ccn_mass),zqi(ig,l,igcm_ccn_number),
464     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
465           
466          ENDDO
467        ENDDO
468       
469       ELSE
470       
471        DO l = 1, nlay
472          DO ig=1,ngrid
473         
474            call updaterice_typ(zqi(ig,l,igcm_h2o_ice),
475     &                      tau(ig,1),zlay(ig,l),rice(ig,l))
476
477          ENDDO
478        ENDDO
479       ENDIF ! of IF(microphys)
480      endif ! of if (water)
481
482      END
483
Note: See TracBrowser for help on using the repository browser.