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

Last change on this file since 3567 was 3360, checked in by emillour, 7 months ago

Generic PCM:
OpenMP fix in newsedim, add a missing THREADPRIVATE for "molrad" variable.
While at it make sedim and callsedim modules.
EM

File size: 5.9 KB
Line 
1      MODULE callsedim_mod
2     
3      IMPLICIT NONE
4     
5      CONTAINS
6     
7      SUBROUTINE callsedim(ngrid,nlay, ptimestep,
8     &                pplev,zlev, pt, pdt,
9     &                pq, pdqfi, pdqsed,pdqs_sed,nq)
10
11      use radinc_h, only : naerkind
12      use radii_mod, only: h2o_reffrad
13      use aerosol_mod, only : iaero_h2o
14      USE tracer_h, only : igcm_co2_ice,igcm_h2o_ice,radius,rho_q
15      use comcstfi_mod, only: g
16      use callkeys_mod, only : water
17      use newsedim_mod, only: newsedim
18
19      IMPLICIT NONE
20
21!==================================================================
22!     
23!     Purpose
24!     -------
25!     Calculates sedimentation of aerosols depending on their
26!     density and radius.
27!     
28!     Authors
29!     -------
30!     F. Forget (1999)
31!     Tracer generalisation by E. Millour (2009)
32!     
33!==================================================================
34
35c-----------------------------------------------------------------------
36c   declarations:
37c   -------------
38
39c   arguments:
40c   ----------
41
42      integer,intent(in):: ngrid ! number of horizontal grid points
43      integer,intent(in):: nlay  ! number of atmospheric layers
44      real,intent(in):: ptimestep  ! physics time step (s)
45      real,intent(in):: pplev(ngrid,nlay+1) ! pressure at inter-layers (Pa)
46      real,intent(in):: pt(ngrid,nlay)      ! temperature at mid-layer (K)
47      real,intent(in):: pdt(ngrid,nlay) ! tendency on temperature
48      real,intent(in):: zlev(ngrid,nlay+1)  ! altitude at layer boundaries
49      integer,intent(in) :: nq ! number of tracers
50      real,intent(in) :: pq(ngrid,nlay,nq)  ! tracers (kg/kg)
51      real,intent(in) :: pdqfi(ngrid,nlay,nq)  ! tendency on tracers before
52                                               ! sedimentation (kg/kg.s-1)
53     
54      real,intent(out) :: pdqsed(ngrid,nlay,nq) ! tendency due to sedimentation (kg/kg.s-1)
55      real,intent(out) :: pdqs_sed(ngrid,nq)    ! flux at surface (kg.m-2.s-1)
56
57c   local:
58c   ------
59
60      INTEGER l,ig, iq
61
62      ! for particles with varying radii:
63      real,allocatable,save :: reffrad(:,:,:) ! particle radius (m)
64      real,allocatable,save :: nueffrad(:,:,:) ! aerosol effective radius variance
65!$OMP THREADPRIVATE(reffrad,nueffrad)
66
67      real zqi(ngrid,nlay,nq) ! to locally store tracers
68      real zt(ngrid,nlay) ! to locally store temperature (K)
69      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
70      real epaisseur (ngrid,nlay) ! Layer thickness (m)
71      real wq(ngrid,nlay+1) ! displaced tracer mass (kg.m-2)
72c      real dens(ngrid,nlay) ! Mean density of the ice part. accounting for dust core
73
74
75      LOGICAL,SAVE :: firstcall=.true.
76!$OMP THREADPRIVATE(firstcall)
77
78c    ** un petit test de coherence
79c       --------------------------
80
81      IF (firstcall) THEN
82        firstcall=.false.
83        ! add some tests on presence of required tracers/aerosols:
84        if (water.and.(igcm_h2o_ice.eq.0)) then
85            write(*,*) "callsedim error: water=.true.",
86     &                 " but igcm_h2o_ice=0"
87          stop
88        endif
89        ! allocate "naerkind" size local arrays (which are also
90        ! "saved" so that this is done only once in for all even if
91        ! we don't need to store the value from a time step to the next)
92        allocate(reffrad(ngrid,nlay,naerkind))
93        allocate(nueffrad(ngrid,nlay,naerkind))
94      ENDIF ! of IF (firstcall)
95     
96!=======================================================================
97!     Preliminary calculation of the layer characteristics
98!     (mass (kg.m-2), thickness (m), etc.
99
100      do  l=1,nlay
101        do ig=1, ngrid
102          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
103          epaisseur(ig,l)= zlev(ig,l+1) - zlev(ig,l)
104          zt(ig,l)=pt(ig,l)+pdt(ig,l)*ptimestep
105        end do
106      end do
107
108!======================================================================
109! Calculate the transport due to sedimentation for each tracer
110! [This has been rearranged by L. Kerber to allow the sedimentation
111!  of general tracers.]
112 
113      do iq=1,nq
114       if((radius(iq).gt.1.e-9).and.(iq.ne.igcm_co2_ice)) then ! JVO 08/2017 : be careful radius was tested uninitialized (fixed) ...
115       
116!         (no sedim for gases, and co2_ice sedim is done in condense_co2)     
117
118! store locally updated tracers
119
120          do l=1,nlay
121            do ig=1, ngrid
122              zqi(ig,l,iq)=pq(ig,l,iq)+pdqfi(ig,l,iq)*ptimestep
123            enddo
124          enddo ! of do l=1,nlay
125       
126!======================================================================
127! Sedimentation
128!======================================================================
129! Water         
130       if (water.and.(iaero_h2o.ne.0).and.(iq.eq.igcm_h2o_ice)) then
131            ! compute radii for h2o_ice
132             call h2o_reffrad(ngrid,nlay,zqi(1,1,igcm_h2o_ice),zt,
133     &                reffrad(1,1,iaero_h2o),nueffrad(1,1,iaero_h2o))
134            ! call sedimentation for h2o_ice
135             call newsedim(ngrid,nlay,ngrid*nlay,ptimestep,
136     &            pplev,masse,epaisseur,zt,reffrad(1,1,iaero_h2o),
137     &            rho_q(iq),zqi(1,1,igcm_h2o_ice),wq,iq)
138
139! General Case
140       else
141             call newsedim(ngrid,nlay,1,ptimestep,
142     &            pplev,masse,epaisseur,zt,radius(iq),rho_q(iq),
143     &            zqi(1,1,iq),wq,iq)
144       endif
145
146!=======================================================================
147!     Calculate the tendencies
148!======================================================================
149
150          do ig=1,ngrid
151!     Ehouarn: with new way of tracking tracers by name, this is simply
152            pdqs_sed(ig,iq) = wq(ig,1)/ptimestep
153          end do
154
155          DO l = 1, nlay
156            DO ig=1,ngrid
157              pdqsed(ig,l,iq)=(zqi(ig,l,iq)-
158     &        (pq(ig,l,iq) + pdqfi(ig,l,iq)*ptimestep))/ptimestep
159            ENDDO
160          ENDDO
161       endif ! of no gases no co2_ice
162      enddo ! of do iq=1,nq
163
164      END SUBROUTINE callsedim
165     
166      END MODULE callsedim_mod
Note: See TracBrowser for help on using the repository browser.