source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/rain.f90 @ 224

Last change on this file since 224 was 224, checked in by ymipsl, 10 years ago
File size: 13.9 KB
Line 
1subroutine rain(ngrid,nq,ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
2
3
4  use ioipsl_getincom, only: getin
5  use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2,Psat_water,Tsat_water,rhowater
6  use radii_mod, only: h2o_cloudrad
7  USE tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
8  implicit none
9
10!==================================================================
11!     
12!     Purpose
13!     -------
14!     Calculates H2O precipitation using simplified microphysics.
15!     
16!     Authors
17!     -------
18!     Adapted from the LMDTERRE code by R. Wordsworth (2009)
19!     Added rain vaporization in case of T>Tsat
20!     Original author Z. X. Li (1993)
21!     
22!==================================================================
23
24  include "dimensions.h"
25  include "dimphys.h"
26  include "comcstfi.h"
27  include "callkeys.h"
28
29!     Arguments
30      integer,intent(in) :: ngrid ! number of atmospherci columns
31      integer,intent(in) :: nq ! number of tracers
32      real,intent(in) :: ptimestep    ! time interval
33      real,intent(in) :: pplev(ngrid,nlayermx+1) ! inter-layer pressure (Pa)
34      real,intent(in) :: pplay(ngrid,nlayermx)   ! mid-layer pressure (Pa)
35      real,intent(in) :: t(ngrid,nlayermx) ! input temperature (K)
36      real,intent(in) :: pdt(ngrid,nlayermx) ! input tendency on temperature (K/s)     
37      real,intent(in) :: pq(ngrid,nlayermx,nq)  ! tracers (kg/kg)
38      real,intent(in) :: pdq(ngrid,nlayermx,nq) ! input tendency on tracers
39      real,intent(out) :: d_t(ngrid,nlayermx) ! temperature tendency (K/s)
40      real,intent(out) :: dqrain(ngrid,nlayermx,nq) ! tendency of H2O precipitation (kg/kg.s-1)
41      real,intent(out) :: dqsrain(ngrid)  ! rain flux at the surface (kg.m-2.s-1)
42      real,intent(out) :: dqssnow(ngrid)  ! snow flux at the surface (kg.m-2.s-1)
43      real,intent(in) :: rneb(ngrid,nlayermx) ! cloud fraction
44
45      REAL zt(ngrid,nlayermx)         ! working temperature (K)
46      REAL ql(ngrid,nlayermx)         ! liquid water (Kg/Kg)
47      REAL q(ngrid,nlayermx)          ! specific humidity (Kg/Kg)
48      REAL d_q(ngrid,nlayermx)        ! water vapor increment
49      REAL d_ql(ngrid,nlayermx)       ! liquid water / ice increment
50
51!     Subroutine options
52      REAL,PARAMETER :: seuil_neb=0.001  ! Nebulosity threshold
53
54      INTEGER,save :: precip_scheme      ! id number for precipitaion scheme
55!     for simple scheme  (precip_scheme=1)
56      REAL,SAVE :: rainthreshold                ! Precipitation threshold in simple scheme
57!     for sundquist scheme  (precip_scheme=2-3)
58      REAL,SAVE :: cloud_sat                    ! Precipitation threshold in non simple scheme
59      REAL,SAVE :: precip_timescale             ! Precipitation timescale
60!     for Boucher scheme  (precip_scheme=4)
61      REAL,SAVE :: Cboucher ! Precipitation constant in Boucher 95 scheme
62      REAL,PARAMETER :: Kboucher=1.19E8
63      REAL,SAVE :: c1
64
65      INTEGER,PARAMETER :: ninter=5
66
67      logical,save :: evap_prec ! Does the rain evaporate?
68
69!     for simple scheme
70      real,parameter :: t_crit=218.0
71      real lconvert
72
73!     Local variables
74      INTEGER i, k, n
75      REAL zqs(ngrid,nlayermx),Tsat(ngrid,nlayermx), zdelta, zcor
76      REAL zrfl(ngrid), zrfln(ngrid), zqev, zqevt
77
78      REAL zoliq(ngrid)
79      REAL zdz(ngrid),zrho(ngrid),ztot(ngrid), zrhol(ngrid)
80      REAL zchau(ngrid),zfroi(ngrid),zfrac(ngrid),zneb(ngrid)
81
82      real reffh2oliq(ngrid,nlayermx),reffh2oice(ngrid,nlayermx)
83 
84      real ttemp, ptemp, psat_tmp
85      real tnext(ngrid,nlayermx)
86
87      real l2c(ngrid,nlayermx)
88      real dWtot
89
90
91!     Indices of water vapour and water ice tracers
92      INTEGER, SAVE :: i_vap=0  ! water vapour
93      INTEGER, SAVE :: i_ice=0  ! water ice
94
95      LOGICAL,SAVE :: firstcall=.true.
96
97!     Online functions
98      REAL fallv, fall2v, zzz ! falling speed of ice crystals
99      fallv (zzz) = 3.29 * ((zzz)**0.16)
100      fall2v (zzz) =10.6 * ((zzz)**0.31)  !for use with radii
101
102
103      IF (firstcall) THEN
104
105         i_vap=igcm_h2o_vap
106         i_ice=igcm_h2o_ice
107       
108         write(*,*) "rain: i_ice=",i_ice
109         write(*,*) "      i_vap=",i_vap
110
111         PRINT*, 'in rain.F, ninter=', ninter
112         PRINT*, 'in rain.F, evap_prec=', evap_prec
113
114         write(*,*) "Precipitation scheme to use?"
115         precip_scheme=1 ! default value
116         call getin("precip_scheme",precip_scheme)
117         write(*,*) " precip_scheme = ",precip_scheme
118
119         if (precip_scheme.eq.1) then
120            write(*,*) "rainthreshold in simple scheme?"
121            rainthreshold=0. ! default value
122            call getin("rainthreshold",rainthreshold)
123            write(*,*) " rainthreshold = ",rainthreshold
124
125         else if (precip_scheme.eq.2.or.precip_scheme.eq.3) then
126            write(*,*) "cloud water saturation level in non simple scheme?"
127            cloud_sat=2.6e-4   ! default value
128            call getin("cloud_sat",cloud_sat)
129            write(*,*) " cloud_sat = ",cloud_sat
130            write(*,*) "precipitation timescale in non simple scheme?"
131            precip_timescale=3600.  ! default value
132            call getin("precip_timescale",precip_timescale)
133            write(*,*) " precip_timescale = ",precip_timescale
134
135         else if (precip_scheme.eq.4) then
136            write(*,*) "multiplicative constant in Boucher 95 precip scheme"
137            Cboucher=1.   ! default value
138            call getin("Cboucher",Cboucher)
139            write(*,*) " Cboucher = ",Cboucher 
140            c1=1.00*1.097/rhowater*Cboucher*Kboucher 
141
142         endif
143
144         write(*,*) "re-evaporate precipitations?"
145         evap_prec=.true. ! default value
146         call getin("evap_prec",evap_prec)
147         write(*,*) " evap_prec = ",evap_prec
148
149         firstcall = .false.
150      ENDIF ! of IF (firstcall)
151
152!     GCM -----> subroutine variables
153      DO k = 1, nlayermx
154      DO i = 1, ngrid
155
156         zt(i,k)   = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here
157         q(i,k)    = pq(i,k,i_vap)+pdq(i,k,i_vap)*ptimestep
158         ql(i,k)   = pq(i,k,i_ice)+pdq(i,k,i_ice)*ptimestep
159
160         !q(i,k)    = pq(i,k,i_vap)!+pdq(i,k,i_vap)
161         !ql(i,k)   = pq(i,k,i_ice)!+pdq(i,k,i_ice)
162
163         if(q(i,k).lt.0.)then ! if this is not done, we don't conserve water
164            q(i,k)=0.
165         endif
166         if(ql(i,k).lt.0.)then
167            ql(i,k)=0.
168         endif
169
170      ENDDO
171      ENDDO
172
173!     Initialise the outputs
174      d_t(1:ngrid,1:nlayermx) = 0.0
175      d_q(1:ngrid,1:nlayermx) = 0.0
176      d_ql(1:ngrid,1:nlayermx) = 0.0
177      zrfl(1:ngrid) = 0.0
178      zrfln(1:ngrid) = 0.0
179
180      ! calculate saturation mixing ratio
181      DO k = 1, nlayermx
182         DO i = 1, ngrid
183            ttemp = zt(i,k) 
184            ptemp = pplay(i,k)
185!            call watersat(ttemp,ptemp,zqs(i,k))
186            call Psat_water(ttemp,ptemp,psat_tmp,zqs(i,k))
187            call Tsat_water(ptemp,Tsat(i,k))
188         ENDDO
189      ENDDO
190
191      ! get column / layer conversion factor
192      DO k = 1, nlayermx
193         DO i = 1, ngrid
194            l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
195         ENDDO
196      ENDDO
197
198      ! Vertical loop (from top to bottom)
199      ! We carry the rain with us and calculate that added by warm/cold precipitation
200      ! processes and that subtracted by evaporation at each level.
201      DO k = nlayermx, 1, -1
202
203         IF (evap_prec) THEN ! note no rneb dependence!
204            DO i = 1, ngrid
205               IF (zrfl(i) .GT.0.) THEN
206
207                  if(zt(i,k).gt.Tsat(i,k))then
208!!                   treat the case where all liquid water should boil
209                     zqev=MIN((zt(i,k)-Tsat(i,k))*RCPD*l2c(i,k)/RLVTT/ptimestep,zrfl(i))
210                     zrfl(i)=MAX(zrfl(i)-zqev,0.)
211                     d_q(i,k)=zqev/l2c(i,k)*ptimestep
212                     d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
213                  else
214                     zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*l2c(i,k)/ptimestep !there was a bug here
215                     zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k))    & !default was 2.e-5
216                        *sqrt(zrfl(i))*l2c(i,k)/pplay(i,k)*zt(i,k)*R ! BC modif here
217                     zqevt = MAX (zqevt, 0.0)
218                     zqev  = MIN (zqev, zqevt)
219                     zqev  = MAX (zqev, 0.0)
220                     zrfln(i)= zrfl(i) - zqev
221                     zrfln(i)= max(zrfln(i),0.0)
222
223                     d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
224                     !d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
225                     d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
226                     zrfl(i)  = zrfln(i)
227                  end if
228                     
229
230               ENDIF ! of IF (zrfl(i) .GT.0.)
231            ENDDO
232         ENDIF ! of IF (evap_prec)
233
234         zoliq(1:ngrid) = 0.0
235
236
237         if(precip_scheme.eq.1)then
238
239            DO i = 1, ngrid
240               ttemp = zt(i,k)
241               IF (ttemp .ge. T_h2O_ice_liq) THEN
242                  lconvert=rainthreshold
243               ELSEIF (ttemp .gt. t_crit) THEN
244                  lconvert=rainthreshold*(1.- t_crit/ttemp)
245                  lconvert=MAX(0.0,lconvert)             
246               ELSE
247                  lconvert=0.
248               ENDIF
249
250
251               IF (ql(i,k).gt.1.e-9) then
252                  zneb(i)  = MAX(rneb(i,k), seuil_neb)
253                  IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
254                     d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
255                     zrfl(i)   = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
256                  ENDIF
257               ENDIF
258            ENDDO
259
260         elseif (precip_scheme.ge.2) then
261         
262           DO i = 1, ngrid
263               IF (rneb(i,k).GT.0.0) THEN
264                  zoliq(i) = ql(i,k)
265                  zrho(i)  = pplay(i,k) / ( zt(i,k) * R )
266                  zdz(i)   = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
267                  zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
268                  zfrac(i) = MAX(zfrac(i), 0.0)
269                  zfrac(i) = MIN(zfrac(i), 1.0)
270                  zneb(i)  = MAX(rneb(i,k), seuil_neb)
271               ENDIF
272           ENDDO
273
274 !recalculate liquid water particle radii
275           call h2o_cloudrad(ngrid,ql,reffh2oliq,reffh2oice)
276
277           SELECT CASE(precip_scheme)
278 !precip scheme from Sundquist 78
279           CASE(2)
280
281            DO n = 1, ninter
282               DO i = 1, ngrid
283                  IF (rneb(i,k).GT.0.0) THEN
284                     ! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
285
286                     zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale)) * zoliq(i)      &
287                          * (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat)**2)) * zfrac(i)
288                     zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
289                     zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
290                          *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly...
291                     ztot(i)  = zchau(i) + zfroi(i)
292
293                     IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
294                     ztot(i)  = MIN(MAX(ztot(i),0.0),zoliq(i))
295                     zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
296
297                  ENDIF
298               ENDDO
299            ENDDO           
300
301 !precip scheme modified from Sundquist 78 (in q**3)
302           CASE(3)         
303           
304            DO n = 1, ninter
305               DO i = 1, ngrid
306                  IF (rneb(i,k).GT.0.0) THEN
307                     ! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
308
309                     zchau(i) = (ptimestep/(FLOAT(ninter)*precip_timescale*cloud_sat**2)) * (zoliq(i)/zneb(i))**3 
310                     zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
311                     zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
312                          *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly...
313                     ztot(i)  = zchau(i) + zfroi(i)
314
315                     IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
316                     ztot(i)  = MIN(MAX(ztot(i),0.0),zoliq(i))
317                     zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
318
319                  ENDIF
320               ENDDO
321            ENDDO           
322
323 !precip scheme modified from Boucher 95
324           CASE(4)
325
326            DO n = 1, ninter
327               DO i = 1, ngrid
328                  IF (rneb(i,k).GT.0.0) THEN
329                     ! this is the ONLY place where zneb and c1 are used
330
331                     zchau(i) = ptimestep/FLOAT(ninter) *c1* zrho(i) &
332                                    *(zoliq(i)/zneb(i))**2*reffh2oliq(i,k)*zneb(i)* zfrac(i) 
333                     zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
334                     zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i)    &
335                          *fall2v(reffh2oice(i,k)) * (1.0-zfrac(i)) ! zfroi behaves oddly...
336                     ztot(i)  = zchau(i) + zfroi(i)
337
338                     IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
339                     ztot(i)  = MIN(MAX(ztot(i),0.0),zoliq(i))
340                     zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
341
342                  ENDIF
343               ENDDO
344            ENDDO           
345
346           END SELECT ! precip_scheme
347
348            ! Change in cloud density and surface H2O values
349            DO i = 1, ngrid
350               IF (rneb(i,k).GT.0.0) THEN
351                  d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
352                  zrfl(i)   = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
353               ENDIF
354            ENDDO
355
356
357         endif ! if precip_scheme=1
358
359      ENDDO ! of DO k = nlayermx, 1, -1
360
361!     Rain or snow on the ground
362      DO i = 1, ngrid
363         if(zrfl(i).lt.0.0)then
364            print*,'Droplets of negative rain are falling...'
365            call abort
366         endif
367         IF (t(i,1) .LT. T_h2O_ice_liq) THEN
368            dqssnow(i) = zrfl(i)
369            dqsrain(i) = 0.0
370         ELSE
371            dqssnow(i) = 0.0
372            dqsrain(i) = zrfl(i) ! liquid water = ice for now
373         ENDIF
374      ENDDO
375
376!     now subroutine -----> GCM variables
377      if (evap_prec) then
378        dqrain(1:ngrid,1:nlayermx,i_vap)=d_q(1:ngrid,1:nlayermx)/ptimestep
379        d_t(1:ngrid,1:nlayermx)=d_t(1:ngrid,1:nlayermx)/ptimestep
380      else
381        dqrain(1:ngrid,1:nlayermx,i_vap)=0.0
382        d_t(1:ngrid,1:nlayermx)=0.0
383      endif
384      dqrain(1:ngrid,1:nlayermx,i_ice) = d_ql(1:ngrid,1:nlayermx)/ptimestep
385
386    end subroutine rain
Note: See TracBrowser for help on using the repository browser.