source: trunk/LMDZ.GENERIC/libf/phystd/callsedim.F @ 1243

Last change on this file since 1243 was 1006, checked in by emillour, 11 years ago

Generic GCM:

  • update the sponge layer: trun it into a module and (more important) compute the sponge quenching analytically rather than via Forward Euler approximation.

EM

File size: 5.5 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#include "dimensions.h"
31#include "dimphys.h"
32#include "comcstfi.h"
33#include "callkeys.h"
34
35c   arguments:
36c   ----------
37
38      integer,intent(in):: ngrid ! number of horizontal grid points
39      integer,intent(in):: nlay  ! number of atmospheric layers
40      real,intent(in):: ptimestep  ! physics time step (s)
41      real,intent(in):: pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
42      real,intent(in):: pt(ngrid,nlay)      ! temperature at mid-layer (K)
43      real,intent(in):: pdt(ngrid,nlay) ! tendency on temperature
44      real,intent(in):: zlev(ngrid,nlay+1)  ! altitude at layer boundaries
45      integer,intent(in) :: nq ! number of tracers
46      real,intent(in) :: pq(ngrid,nlay,nq)  ! tracers (kg/kg)
47      real,intent(in) :: pdqfi(ngrid,nlay,nq)  ! tendency on tracers before
48                                               ! sedimentation (kg/kg.s-1)
49     
50      real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
51      real,intent(out) :: pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)
52     
53c   local:
54c   ------
55
56      INTEGER l,ig, iq
57
58      ! for particles with varying radii:
59      real reffrad(ngrid,nlay,naerkind) ! particle radius (m)
60      real nueffrad(ngrid,nlay,naerkind) ! aerosol effective radius variance
61
62      real zqi(ngrid,nlay,nq) ! to locally store tracers
63      real zt(ngrid,nlay) ! to locally store temperature (K)
64      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
65      real epaisseur (ngrid,nlay) ! Layer thickness (m)
66      real wq(ngrid,nlay+1) ! displaced tracer mass (kg.m-2)
67c      real dens(ngrid,nlayermx) ! Mean density of the ice part. accounting for dust core
68
69
70      LOGICAL,SAVE :: firstcall=.true.
71
72c    ** un petit test de coherence
73c       --------------------------
74
75      IF (firstcall) THEN
76        firstcall=.false.
77        ! add some tests on presence of required tracers/aerosols:
78        if (water) then
79          if (igcm_h2o_ice.eq.0) then
80            write(*,*) "callsedim error: water=.true.",
81     &                 " but igcm_h2o_ice=0"
82          stop
83          endif
84          if (iaero_h2o.eq.0) then
85            write(*,*) "callsedim error: water=.true.",
86     &                 " but iaero_ho2=0"
87          stop
88          endif
89        endif
90      ENDIF ! of IF (firstcall)
91     
92!=======================================================================
93!     Preliminary calculation of the layer characteristics
94!     (mass (kg.m-2), thickness (m), etc.
95
96      do  l=1,nlay
97        do ig=1, ngrid
98          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
99          epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
100          zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
101        end do
102      end do
103
104!======================================================================
105! Calculate the transport due to sedimentation for each tracer
106! [This has been rearranged by L. Kerber to allow the sedimentation
107!  of general tracers.]
108 
109      do iq=1,nq
110       if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_co2_ice)) then   
111!         (no sedim for gases, and co2_ice sedim is done in condense_cloud)     
112
113! store locally updated tracers
114
115          do l=1,nlay
116            do ig=1, ngrid
117              zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
118            enddo
119          enddo ! of do l=1,nlay
120       
121!======================================================================
122! Sedimentation
123!======================================================================
124! Water         
125          if (water.and.(iq.eq.igcm_h2o_ice)) then
126            ! compute radii for h2o_ice
127             call h2o_reffrad(ngrid,zqi(1,1,igcm_h2o_ice),zt,
128     &                reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
129            ! call sedimentation for h2o_ice
130             call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
131     &            pplev,masse,epaisseur,zt,reffrad(1,1,iaero_h2o),
132     &            rho_q(iq),zqi(1,1,igcm_h2o_ice),wq)
133
134! General Case
135          else
136             call newsedim(ngrid,nlay,1,ptimestep,
137     &            pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
138     &            zqi(1,1,iq),wq)
139          endif
140
141!=======================================================================
142!     Calculate the tendencies
143!======================================================================
144
145          do ig=1,ngrid
146!     Ehouarn: with new way of tracking tracers by name, this is simply
147            pdqs_sed(ig,iq) = wq(ig,1)/ptimestep
148          end do
149
150          DO l = 1, nlay
151            DO ig=1,ngrid
152              pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
153     &        (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
154            ENDDO
155          ENDDO
156       endif ! of no gases no co2_ice
157      enddo ! of do iq=1,nq
158      return
159      end
Note: See TracBrowser for help on using the repository browser.