source: trunk/WRF.COMMON/WRFV2/phys/module_mp_kessler.F @ 3094

Last change on this file since 3094 was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien

File size: 8.9 KB
Line 
1!WRF:MODEL_LAYER:PHYSICS
2!
3
4MODULE module_mp_kessler
5
6CONTAINS
7!----------------------------------------------------------------
8   SUBROUTINE kessler( t, qv, qc, qr, rho, pii                  &
9                      ,dt_in, z, xlv, cp                        &
10                      ,EP2,SVP1,SVP2,SVP3,SVPT0,rhowater        &
11                      ,dz8w                                     &
12                      ,RAINNC, RAINNCV                          &
13                      ,ids,ide, jds,jde, kds,kde                & ! domain dims
14                      ,ims,ime, jms,jme, kms,kme                & ! memory dims
15                      ,its,ite, jts,jte, kts,kte                & ! tile   dims
16                                                                )
17!----------------------------------------------------------------
18   IMPLICIT NONE
19!----------------------------------------------------------------
20   !  taken from the COMMAS code - WCS 10 May 1999.
21   !  converted from FORTRAN 77 to 90, tiled, WCS 10 May 1999.
22!----------------------------------------------------------------
23   REAL    , PARAMETER ::  c1 = .001
24   REAL    , PARAMETER ::  c2 = .001
25   REAL    , PARAMETER ::  c3 = 2.2
26   REAL    , PARAMETER ::  c4 = .875
27   REAL    , PARAMETER ::  fudge = 1.0
28   REAL    , PARAMETER ::  mxfall = 10.0
29!----------------------------------------------------------------
30   INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde, &
31                                     ims,ime, jms,jme, kms,kme, &
32                                     its,ite, jts,jte, kts,kte
33   REAL   ,      INTENT(IN   )    :: xlv, cp
34   REAL   ,      INTENT(IN   )    :: EP2,SVP1,SVP2,SVP3,SVPT0
35   REAL   ,      INTENT(IN   )    :: rhowater
36
37   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),              &
38         INTENT(INOUT) ::                                       &
39                                                            t , &
40                                                            qv, &
41                                                            qc, &
42                                                            qr
43
44   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),              &
45         INTENT(IN   ) ::                                       &
46                                                           rho, &
47                                                           pii, &
48                                                          dz8w
49
50   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ),              &
51         INTENT(IN   ) ::                                    z
52
53   REAL, INTENT(IN   ) :: dt_in
54
55   REAL, DIMENSION( ims:ime , jms:jme ),                        &
56         INTENT(INOUT) ::                               RAINNC, &
57                                                       RAINNCV
58
59
60
61   ! local variables
62
63   REAL :: qrprod, ern, gam, rcgs, rcgsi
64   REAL, DIMENSION( its:ite , kts:kte, jts:jte ) ::     prod
65   REAL, DIMENSION(kts:kte) :: vt, prodk, vtden,rdzk,rhok,factor,rdzw
66   INTEGER :: i,j,k
67   INTEGER :: nfall, n, nfall_new
68   REAL    :: qrr, pressure, temp, es, qvs, dz, dt
69   REAL    :: f5, dtfall, rdz, product
70   REAL    :: max_heating, max_condense, max_rain, maxqrp
71   REAL    :: vtmax, ernmax, crmax, factorn, time_sediment
72   REAL    :: qcr, factorr, ppt
73   REAL, PARAMETER :: max_cr_sedimentation = 0.75
74!----------------------------------------------------------------
75
76   INTEGER :: imax, kmax
77
78    dt = dt_in
79
80!   f5 = 237.3 * 17.27 * 2.5e6 / cp
81    f5 = svp2*(svpt0-svp3)*xlv/cp
82    ernmax = 0.
83    maxqrp = -100.
84
85!------------------------------------------------------------------------------
86! parameters for the time split terminal advection
87!------------------------------------------------------------------------------
88
89      max_heating = 0.
90      max_condense = 0.
91      max_rain = 0.
92
93!-----------------------------------------------------------------------------
94! outer J loop for entire microphysics, outer i loop for sedimentation
95!-----------------------------------------------------------------------------
96
97  microphysics_outer_j_loop: DO j = jts, jte
98
99  sedimentation_outer_i_loop: DO i = its,ite
100
101!   vtmax = 0.
102   crmax = 0.
103
104
105!------------------------------------------------------------------------------
106! Terminal velocity calculation and advection, set up coefficients and
107! compute stable timestep
108!------------------------------------------------------------------------------
109
110   DO k = 1, kte
111     prodk(k)   = qr(i,k,j)
112     rhok(k) = rho(i,k,j)
113     qrr = amax1(0.,qr(i,k,j)*0.001*rhok(k))
114     vtden(k) = sqrt(rhok(1)/rhok(k))
115     vt(k) = 36.34*(qrr**0.1364) * vtden(k)
116!     vtmax = amax1(vt(k), vtmax)
117     rdzw(k) = 1./dz8w(i,k,j)
118     crmax = amax1(vt(k)*dt*rdzw(k),crmax)
119   ENDDO
120   DO k = 1, kte-1
121     rdzk(k) = 1./(z(i,k+1,j) - z(i,k,j))
122   ENDDO
123   rdzk(kte) = 1./(z(i,kte,j) - z(i,kte-1,j))
124
125   nfall = max(1,nint(0.5+crmax/max_cr_sedimentation))  ! courant number for big timestep.
126   dtfall = dt / float(nfall)                           ! splitting so courant number for sedimentation
127   time_sediment = dt                                   ! is stable
128
129!------------------------------------------------------------------------------
130! Terminal velocity calculation and advection
131! Do a time split loop on this for stability.
132!------------------------------------------------------------------------------
133
134   column_sedimentation: DO WHILE ( nfall > 0 )
135
136   time_sediment = time_sediment - dtfall
137   DO k = 1, kte-1
138     factor(k) = dtfall*rdzk(k)/rhok(k)
139   ENDDO
140   factor(kte) = dtfall*rdzk(kte)
141
142   ppt=0.
143
144      k = 1
145      ppt=rhok(k)*prodk(k)*vt(k)*dtfall/rhowater
146      RAINNCV(i,j)=ppt*1000.
147      RAINNC(i,j)=RAINNC(i,j)+ppt*1000.  ! unit = mm
148 
149!------------------------------------------------------------------------------
150! Time split loop, Fallout done with flux upstream
151!------------------------------------------------------------------------------
152
153      DO k = kts, kte-1
154        prodk(k) = prodk(k) - factor(k)           &
155                  * (rhok(k)*prodk(k)*vt(k)       &
156                    -rhok(k+1)*prodk(k+1)*vt(k+1))
157      ENDDO
158
159      k = kte
160      prodk(k) = prodk(k) - factor(k)*prodk(k)*vt(k)
161
162!------------------------------------------------------------------------------
163! compute new sedimentation velocity, and check/recompute new
164! sedimentation timestep if this isn't the last split step.
165!------------------------------------------------------------------------------
166
167      IF( nfall > 1 ) THEN ! this wasn't the last split sedimentation timestep
168
169        nfall = nfall - 1
170        crmax = 0.
171        DO k = kts, kte
172          qrr = amax1(0.,prodk(k)*0.001*rhok(k))
173          vt(k) = 36.34*(qrr**0.1364) * vtden(k)
174!          vtmax = amax1(vt(k), vtmax)
175          crmax = amax1(vt(k)*time_sediment*rdzw(k),crmax)
176        ENDDO
177
178        nfall_new = max(1,nint(0.5+crmax/max_cr_sedimentation))
179        if (nfall_new /= nfall ) then
180          nfall = nfall_new
181          dtfall = time_sediment/nfall
182        end if
183
184      ELSE  ! this was the last timestep
185
186        DO k=kts,kte
187          prod(i,k,j) = prodk(k)
188        ENDDO
189        nfall = 0  ! exit condition for sedimentation loop
190
191      END IF
192
193   ENDDO column_sedimentation
194
195   ENDDO sedimentation_outer_i_loop
196
197!------------------------------------------------------------------------------
198! Production of rain and deletion of qc
199! Production of qc from supersaturation
200! Evaporation of QR
201!------------------------------------------------------------------------------
202
203     DO k = kts, kte
204     DO i = its, ite
205       factorn = 1.0 / (1.+c3*dt*amax1(0.,qr(i,k,j))**c4)
206       qrprod = qc(i,k,j) * (1.0 - factorn)           &
207             + factorn*c1*dt*amax1(qc(i,k,j)-c2,0.)     
208       rcgs = 0.001*rho(i,k,j)
209
210       qc(i,k,j) = amax1(qc(i,k,j) - qrprod,0.)
211       qr(i,k,j) = (qr(i,k,j) + prod(i,k,j)-qr(i,k,j))
212       qr(i,k,j) = amax1(qr(i,k,j) + qrprod,0.)
213
214       temp      = pii(i,k,j)*t(i,k,j)
215       pressure = 1.000e+05 * (pii(i,k,j)**(1004./287.))
216       gam = 2.5e+06/(1004.*pii(i,k,j))
217!      qvs       = 380.*exp(17.27*(temp-273.)/(temp- 36.))/pressure
218       es        = 1000.*svp1*exp(svp2*(temp-svpt0)/(temp-svp3))
219       qvs       = ep2*es/(pressure-es)
220!      prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+qvs*f5/(temp-36.)**2)
221       prod(i,k,j) = (qv(i,k,j)-qvs) / (1.+pressure/(pressure-es)*qvs*f5/(temp-svp3)**2)
222       ern  = amin1(dt*(((1.6+124.9*(rcgs*qr(i,k,j))**.2046)   &
223          *(rcgs*qr(i,k,j))**.525)/(2.55e8/(pressure*qvs)       &
224          +5.4e5))*(dim(qvs,qv(i,k,j))/(rcgs*qvs)),             &
225          amax1(-prod(i,k,j)-qc(i,k,j),0.),qr(i,k,j))
226
227! Update all variables
228
229       product = amax1(prod(i,k,j),-qc(i,k,j))
230       t (i,k,j) = t(i,k,j) + gam*(product - ern)
231       qv(i,k,j) = amax1(qv(i,k,j) - product + ern,0.)
232       qc(i,k,j) =       qc(i,k,j) + product
233       qr(i,k,j) = qr(i,k,j) - ern
234
235     ENDDO
236     ENDDO
237
238  ENDDO  microphysics_outer_j_loop
239
240  RETURN
241
242  END SUBROUTINE kessler
243
244END MODULE module_mp_kessler
Note: See TracBrowser for help on using the repository browser.