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

Last change on this file since 1330 was 1315, checked in by milmd, 11 years ago

LMDZ.GENERIC. OpenMP directives added in generic physic. When running in pure OpenMP or hybrid OpenMP/MPI, may have some bugs with condense_cloud and wstats routines.

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,nlay) ! Mean density of the ice part. accounting for dust core
68
69
70      LOGICAL,SAVE :: firstcall=.true.
71!$OMP THREADPRIVATE(firstcall)
72
73c    ** un petit test de coherence
74c       --------------------------
75
76      IF (firstcall) THEN
77        firstcall=.false.
78        ! add some tests on presence of required tracers/aerosols:
79        if (water) then
80          if (igcm_h2o_ice.eq.0) then
81            write(*,*) "callsedim error: water=.true.",
82     &                 " but igcm_h2o_ice=0"
83          stop
84          endif
85          if (iaero_h2o.eq.0) then
86            write(*,*) "callsedim error: water=.true.",
87     &                 " but iaero_ho2=0"
88          stop
89          endif
90        endif
91      ENDIF ! of IF (firstcall)
92     
93!=======================================================================
94!     Preliminary calculation of the layer characteristics
95!     (mass (kg.m-2), thickness (m), etc.
96
97      do  l=1,nlay
98        do ig=1, ngrid
99          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
100          epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
101          zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
102        end do
103      end do
104
105!======================================================================
106! Calculate the transport due to sedimentation for each tracer
107! [This has been rearranged by L. Kerber to allow the sedimentation
108!  of general tracers.]
109 
110      do iq=1,nq
111       if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_co2_ice)) then   
112!         (no sedim for gases, and co2_ice sedim is done in condense_cloud)     
113
114! store locally updated tracers
115
116          do l=1,nlay
117            do ig=1, ngrid
118              zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
119            enddo
120          enddo ! of do l=1,nlay
121       
122!======================================================================
123! Sedimentation
124!======================================================================
125! Water         
126          if (water.and.(iq.eq.igcm_h2o_ice)) then
127            ! compute radii for h2o_ice
128             call h2o_reffrad(ngrid,nlay,zqi(1,1,igcm_h2o_ice),zt,
129     &                reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
130            ! call sedimentation for h2o_ice
131             call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
132     &            pplev,masse,epaisseur,zt,reffrad(1,1,iaero_h2o),
133     &            rho_q(iq),zqi(1,1,igcm_h2o_ice),wq)
134
135! General Case
136          else
137             call newsedim(ngrid,nlay,1,ptimestep,
138     &            pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
139     &            zqi(1,1,iq),wq)
140          endif
141
142!=======================================================================
143!     Calculate the tendencies
144!======================================================================
145
146          do ig=1,ngrid
147!     Ehouarn: with new way of tracking tracers by name, this is simply
148            pdqs_sed(ig,iq) = wq(ig,1)/ptimestep
149          end do
150
151          DO l = 1, nlay
152            DO ig=1,ngrid
153              pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
154     &        (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
155            ENDDO
156          ENDDO
157       endif ! of no gases no co2_ice
158      enddo ! of do iq=1,nq
159      return
160      end
Note: See TracBrowser for help on using the repository browser.