source: trunk/LMDZ.MARS/libf/aeronomars/surfacearea.F @ 1038

Last change on this file since 1038 was 1036, checked in by emillour, 12 years ago

Mars GCM: (a first step towards using parallel dynamics)

  • IMPORTANT CHANGE: Implemented dynamic tracers. It is no longer necessary to compile the model with the '-t #' option, number of tracers is simply read from tracer.def file (as before). Adapted makegcm_* scripts (and co.) accordingly. Technical aspects of the switch to dynamic tracers are:
    • advtrac.h (in dyn3d) removed and replaced by module infotrac.F
    • tracer.h (in phymars) removed and replaced by module tracer_mod.F90 (which contains nqmx, the number of tracers, etc. and can be used anywhere in the physics).
  • Included some side cleanups: removed unused files (in dyn3d) anldoppler2.F, anl_mcdstats.F and anl_stats-diag.F, and all the unecessary dimensions.* files in grid/dimension.
  • Checked that changes are clean and that GCM yields identical results (in debug mode) to previous svn version.

EM
File size: 5.4 KB
Line 
1      subroutine surfacearea(ngrid, nlay, ptimestep,
2     $                       pplay, pzlay,
3     $                       pt, pq, pdq, nq, rdust, rice, tau,
4     $                       tauscaling,
5     $                       surfdust, surfice)
6
7      use tracer_mod, only: nuice_sed, igcm_dust_number,
8     &                      igcm_ccn_number, varian, ccn_factor
9      implicit none
10
11!==========================================================================
12!     calculation of the ice and dust surface area (m2/m3)
13!     available for heterogeneous reactions
14!
15!     Franck Lefevre
16!     version 1.2 april 2012
17!==========================================================================
18
19#include "dimensions.h"
20#include "dimphys.h"
21#include "comcstfi.h"
22#include "callkeys.h"
23!#include "tracer.h"
24#include "dimradmars.h"
25#include "chimiedata.h"
26#include "conc.h"
27
28! input
29
30      integer,intent(in) :: ngrid, nlay
31      integer,intent(in) :: nq               ! number of tracers
32      real,intent(in) :: ptimestep           ! physics time step (s)
33      real,intent(in) :: pplay(ngrid,nlay)   ! pressure at mid-layers (Pa)
34      real,intent(in) :: pzlay(ngrid,nlay)   ! altitude at mid-layers (m)
35      real,intent(in) :: pt(ngrid,nlay)      ! temperature at mid-layers (K)
36      real,intent(in) :: pq(ngrid,nlay,nq)   ! tracers (kg/kg)
37      real,intent(in) :: pdq(ngrid,nlay,nq)  ! physical tendency (kg/kg.s-1)
38      real,intent(in) :: rdust(ngrid,nlay)   ! dust geometric mean radius (m)
39      real,intent(in) :: rice(ngrid,nlay)    ! ice mass mean radius (m)
40      real,intent(in) :: tau(ngrid,naerkind) ! column dust optical depth at each point
41      real,intent(in) :: tauscaling(ngrid)   ! conversion factor for dust amount
42
43! output
44
45      real,intent(out) :: surfdust(ngrid,nlay) ! dust surface area (m2/m3)
46      real,intent(out) :: surfice(ngrid,nlay)  ! water-ice surface area (m2/m3)
47
48! local
49
50      integer    :: l, ig
51      real       :: rho                     ! density (kg/m3)
52      real       :: dustnd, icend           ! uodated dust and ice number densities (kg/kg)
53      real, save :: factor_ice, factor_dust ! multiplying factor to compute total surface area
54                                            ! from the mass-mean radius
55      real       :: sigma_ice, sigma_dust   ! variance of the ice and dust distributions
56      real       :: ccntyp                  ! typical dust number density (#/kg)
57                                            ! (microphys = false)
58      real       :: rdusttyp                ! typical dust radius (m)
59                                            ! (microphys = false)
60
61      logical, save :: firstcall = .true.
62
63!==========================================================================
64
65      if (firstcall) then ! compute the multiplying factors
66         sigma_dust  = varian
67         sigma_ice   = sqrt(log(nuice_sed + 1.))
68         factor_dust = exp(0.5*(log(sigma_dust))**2)
69         factor_ice  = exp(0.5*(log(sigma_ice))**2)
70         write(*,*) 'surfacearea : factor_dust = ', factor_dust
71         write(*,*) 'surfacearea : factor_ice  = ', factor_ice
72         firstcall = .false.
73      end if
74
75      if (microphys) then ! improvedclouds
76         do l = 1,nlay
77            do ig = 1,ngrid
78!              atmospheric density
79               rho = pplay(ig,l)/(rnew(ig,l)*pt(ig,l))
80!              updated dust number density
81               dustnd = pq(ig,l,igcm_dust_number)
82     $                + pdq(ig,l,igcm_dust_number)*ptimestep
83!              updated ice number density
84               icend  = pq(ig,l,igcm_ccn_number)
85     $                + pdq(ig,l,igcm_ccn_number)*ptimestep
86!              dust surface area
87               surfdust(ig,l) = factor_dust*dustnd*rho*tauscaling(ig)
88     $                          *4.*pi*rdust(ig,l)**2
89!              ice surface area
90               surfice(ig,l)  = factor_ice*icend*rho*tauscaling(ig)
91     $                          *4.*pi*rice(ig,l)**2
92            end do
93         end do
94      else               ! simpleclouds
95         do l = 1,nlay
96            do ig = 1,ngrid
97!              atmospheric density
98               rho = pplay(ig,l)/(rnew(ig,l)*pt(ig,l))
99!              typical dust radius
100               rdusttyp = max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9)
101!              typical dust number density
102               ccntyp = 1.3e+8*max(tau(ig,1),0.001)/0.1
103     $                  *exp(-pzlay(ig,l)/10000.)
104               ccntyp = ccntyp/ccn_factor
105               if (rice(ig,l) .gt. rdust(ig,l)) then
106                  surfdust(ig,l) = factor_dust*ccntyp*(ccn_factor - 1.)
107     $                             *rho*4.*pi*rdusttyp**2
108                  surfice(ig,l)  = factor_ice*ccntyp*4.*pi*rice(ig,l)**2
109               else
110                  surfdust(ig,l) = factor_dust*ccntyp*ccn_factor
111     $                             *rho*4.*pi*rdusttyp**2
112                  surfice(ig,l)  = 0.
113               end if
114            end do
115         end do
116      end if         ! of microphys
117
118! write diagnostics in micron2/cm3
119     
120      if (callstats) then
121        call wstats(ngrid,"surfdust", "Dust surface area",
122     $            "micron2 cm-3",3,surfdust*1.e6)
123        call wstats(ngrid,"surfice", "Ice cloud surface area",
124     $            "micron2 cm-3",3,surfice*1.e6)
125      endif
126      call writediagfi(ngrid,"surfdust", "Dust surface area",
127     $            "micron2 cm-3",3,surfdust*1.e6)
128      call writediagfi(ngrid,"surfice", "Ice cloud surface area",
129     $            "micron2 cm-3",3,surfice*1.e6)
130
131      return
132      end
Note: See TracBrowser for help on using the repository browser.