source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/callsedim.f @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 11.4 KB
Line 
1      SUBROUTINE callsedim(ngrid,nlay, ptimestep,
2     &                pplev,zlev, pt, pdt,
3     &                pq, pdqfi, pdqsed,pdqs_sed,nq)
4
5      use radinc_h, only : naerkind
6      use radii_mod, only: h2o_reffrad
7      use aerosol_mod, only : iaero_h2o
8      USE tracer_h, only : igcm_co2_ice,igcm_h2o_ice,radius,rho_q
9
10      IMPLICIT NONE
11
12!==================================================================
13!     
14!     Purpose
15!     -------
16!     Calculates sedimentation of aerosols depending on their
17!     density and radius.
18!     
19!     Authors
20!     -------
21!     F. Forget (1999)
22!     Tracer generalisation by E. Millour (2009)
23!     
24!==================================================================
25
26c-----------------------------------------------------------------------
27c   declarations:
28c   -------------
29
30!-----------------------------------------------------------------------
31!   INCLUDE 'dimensions.h'
32!
33!   dimensions.h contient les dimensions du modele
34!   ndm est tel que iim=2**ndm
35!-----------------------------------------------------------------------
36
37      INTEGER iim,jjm,llm,ndm
38
39      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)
40
41!-----------------------------------------------------------------------
42!-----------------------------------------------------------------------
43!   INCLUDE 'dimphys.h'
44
45! ngridmx : number of horizontal grid points
46! note: the -1/jjm term will be 0; unless jj=1
47      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)   
48! nlayermx : number of atmospheric layers
49      integer, parameter :: nlayermx = llm 
50! nsoilmx : number of subterranean layers ! nsoilmx is now in comsoil_h
51      !integer, parameter :: nsoilmx = 4 ! for a test
52      !integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case
53      !integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
54!-----------------------------------------------------------------------
55
56!-----------------------------------------------------------------------
57! INCLUDE "comcstfi.h"
58
59      common/comcstfi/pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
60      common/comcstfi/avocado!,molrad,visc
61     
62      real pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
63      real avocado!,molrad,visc
64
65!
66! For Fortran 77/Fortran 90 compliance always use line continuation
67! symbols '&' in columns 73 and 6
68!
69! Group commons according to their type for minimal performance impact
70
71      COMMON/callkeys_l/callrad,corrk,calldifv,UseTurbDiff,calladj      &
72     &   , co2cond,callsoil                                             &
73     &   , season,diurnal,tlocked,rings_shadow,lwrite                   &
74     &   , callstats,calleofdump                                        &
75     &   , enertest                                                     &
76     &   , callgasvis,continuum,H2Ocont_simple,graybody                 &
77     &   , radfixed                                                     &
78     &   , meanOLR, specOLR                                             &
79     &   , kastprof                                                     &
80     &   , nosurf, oblate                                               &     
81     &   , newtonian, testradtimes                                      &
82     &   , check_cpp_match, force_cpp                                   &
83     &   , rayleigh                                                     &
84     &   , stelbbody                                                    &
85     &   , nearco2cond                                                  &
86     &   , tracer, mass_redistrib, varactive, varfixed                  &
87     &   , sedimentation,water,watercond,waterrain                      &
88     &   , aeroco2,aeroh2o,aeroh2so4,aeroback2lay                       &
89     &   , aerofixco2,aerofixh2o                                        &
90     &   , hydrology, sourceevol                                        &
91     &   , CLFvarying                                                   &
92     &   , strictboundcorrk                                             &                                       
93     &   , ok_slab_ocean                                                &
94     &   , ok_slab_sic                                                  &
95     &   , ok_slab_heat_transp                                         
96
97
98      COMMON/callkeys_i/iaervar,iddist,iradia,startype
99     
100      COMMON/callkeys_r/topdustref,Nmix_co2,dusttau,Fat1AU,stelTbb,     &
101     &                  Tstrat,tplanet,obs_tau_col_tropo,               &
102     &                  obs_tau_col_strato,pres_bottom_tropo,           &
103     &                  pres_top_tropo,pres_bottom_strato,              &
104     &                  pres_top_strato,size_tropo,size_strato,satval,  &
105     &                  CLFfixval,n2mixratio,co2supsat,pceil,albedosnow,&
106     &                  maxicethick,Tsaldiff,tau_relax,cloudlvl,        &
107     &                  icetstep,intheat,flatten,Rmean,J2,MassPlanet
108     
109      logical callrad,corrk,calldifv,UseTurbDiff                        &
110     &   , calladj,co2cond,callsoil                                     &
111     &   , season,diurnal,tlocked,rings_shadow,lwrite                   &
112     &   , callstats,calleofdump                                        &
113     &   , callgasvis,continuum,H2Ocont_simple,graybody                 &
114     &   , strictboundcorrk                                             
115
116      logical enertest
117      logical nonideal
118      logical meanOLR
119      logical specOLR
120      logical kastprof
121      logical newtonian
122      logical check_cpp_match
123      logical force_cpp
124      logical testradtimes
125      logical rayleigh
126      logical stelbbody
127      logical ozone
128      logical nearco2cond
129      logical tracer
130      logical mass_redistrib
131      logical varactive
132      logical varfixed
133      logical radfixed
134      logical sedimentation
135      logical water,watercond,waterrain
136      logical aeroco2,aeroh2o,aeroh2so4,aeroback2lay
137      logical aerofixco2,aerofixh2o
138      logical hydrology
139      logical sourceevol
140      logical CLFvarying
141      logical nosurf
142      logical oblate
143      logical ok_slab_ocean
144      logical ok_slab_sic
145      logical ok_slab_heat_transp
146
147      integer iddist
148      integer iaervar
149      integer iradia
150      integer startype
151
152      real topdustref
153      real Nmix_co2
154      real dusttau
155      real Fat1AU
156      real stelTbb
157      real Tstrat
158      real tplanet
159      real obs_tau_col_tropo
160      real obs_tau_col_strato
161      real pres_bottom_tropo
162      real pres_top_tropo
163      real pres_bottom_strato
164      real pres_top_strato
165      real size_tropo
166      real size_strato
167      real satval
168      real CLFfixval
169      real n2mixratio
170      real co2supsat
171      real pceil
172      real albedosnow
173      real maxicethick
174      real Tsaldiff
175      real tau_relax
176      real cloudlvl
177      real icetstep
178      real intheat
179      real flatten
180      real Rmean
181      real J2
182      real MassPlanet
183
184c   arguments:
185c   ----------
186
187      integer,intent(in):: ngrid ! number of horizontal grid points
188      integer,intent(in):: nlay  ! number of atmospheric layers
189      real,intent(in):: ptimestep  ! physics time step (s)
190      real,intent(in):: pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
191      real,intent(in):: pt(ngrid,nlay)      ! temperature at mid-layer (K)
192      real,intent(in):: pdt(ngrid,nlay) ! tendency on temperature
193      real,intent(in):: zlev(ngrid,nlay+1)  ! altitude at layer boundaries
194      integer,intent(in) :: nq ! number of tracers
195      real,intent(in) :: pq(ngrid,nlay,nq)  ! tracers (kg/kg)
196      real,intent(in) :: pdqfi(ngrid,nlay,nq)  ! tendency on tracers before
197                                               ! sedimentation (kg/kg.s-1)
198     
199      real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
200      real,intent(out) :: pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)
201     
202c   local:
203c   ------
204
205      INTEGER l,ig, iq
206
207      ! for particles with varying radii:
208      real reffrad(ngrid,nlay,naerkind) ! particle radius (m)
209      real nueffrad(ngrid,nlay,naerkind) ! aerosol effective radius variance
210
211      real zqi(ngrid,nlay,nq) ! to locally store tracers
212      real zt(ngrid,nlay) ! to locally store temperature (K)
213      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
214      real epaisseur (ngrid,nlay) ! Layer thickness (m)
215      real wq(ngrid,nlay+1) ! displaced tracer mass (kg.m-2)
216c      real dens(ngrid,nlayermx) ! Mean density of the ice part. accounting for dust core
217
218
219      LOGICAL,SAVE :: firstcall=.true.
220
221c    ** un petit test de coherence
222c       --------------------------
223
224      IF (firstcall) THEN
225        firstcall=.false.
226        ! add some tests on presence of required tracers/aerosols:
227        if (water) then
228          if (igcm_h2o_ice.eq.0) then
229            write(*,*) "callsedim error: water=.true.",
230     &                 " but igcm_h2o_ice=0"
231          stop
232          endif
233          if (iaero_h2o.eq.0) then
234            write(*,*) "callsedim error: water=.true.",
235     &                 " but iaero_ho2=0"
236          stop
237          endif
238        endif
239      ENDIF ! of IF (firstcall)
240     
241!=======================================================================
242!     Preliminary calculation of the layer characteristics
243!     (mass (kg.m-2), thickness (m), etc.
244
245      do  l=1,nlay
246        do ig=1, ngrid
247          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 
248          epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
249          zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
250        end do
251      end do
252
253!======================================================================
254! Calculate the transport due to sedimentation for each tracer
255! [This has been rearranged by L. Kerber to allow the sedimentation
256!  of general tracers.]
257 
258      do iq=1,nq
259       if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_co2_ice)) then   
260!         (no sedim for gases, and co2_ice sedim is done in condense_cloud)     
261
262! store locally updated tracers
263
264          do l=1,nlay 
265            do ig=1, ngrid
266              zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
267            enddo
268          enddo ! of do l=1,nlay
269       
270!======================================================================
271! Sedimentation
272!======================================================================
273! Water         
274          if (water.and.(iq.eq.igcm_h2o_ice)) then
275            ! compute radii for h2o_ice
276             call h2o_reffrad(ngrid,zqi(1,1,igcm_h2o_ice),zt,
277     &                reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
278            ! call sedimentation for h2o_ice
279             call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
280     &            pplev,masse,epaisseur,zt,reffrad(1,1,iaero_h2o),
281     &            rho_q(iq),zqi(1,1,igcm_h2o_ice),wq)
282
283! General Case
284          else
285             call newsedim(ngrid,nlay,1,ptimestep,
286     &            pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
287     &            zqi(1,1,iq),wq)
288          endif
289
290!=======================================================================
291!     Calculate the tendencies
292!======================================================================
293
294          do ig=1,ngrid
295!     Ehouarn: with new way of tracking tracers by name, this is simply
296            pdqs_sed(ig,iq) = wq(ig,1)/ptimestep
297          end do
298
299          DO l = 1, nlay
300            DO ig=1,ngrid
301              pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
302     &        (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
303            ENDDO
304          ENDDO
305       endif ! of no gases no co2_ice
306      enddo ! of do iq=1,nq
307      return
308      end
Note: See TracBrowser for help on using the repository browser.