source: trunk/LMDZ.GENERIC/libf/phystd/rain.F90 @ 848

Last change on this file since 848 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

  1. ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
  2. if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

  • Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
  • RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
  • RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
  • RUN RCM1D : perfect match in normal mode.
  • RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

  • Additional arguments to some subroutines (ngrid and nq)
  • F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

  • F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
  • Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
  • Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
  • Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

  • Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

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