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

Last change on this file since 1930 was 1913, checked in by emillour, 8 years ago

Mars GCM:

  • Forgotten in previous commit: gwprofil.F -> gwprofil_mod.F (here also the

size of an argument, rho, was incorrect in caller orodrag).

  • Turned newsedim.F into a module newsedim_mod.F
  • Adapted co2cloud.F and improvedCO2clouds.F to not use "newunit" to open file

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