source: trunk/LMDZ.MARS/libf/phymars/watercloud.F @ 706

Last change on this file since 706 was 703, checked in by tnavarro, 12 years ago

corrected bug of tracers handling in watercloud that led to a incorrect CO2 depletion. Actually, all non-water cycle tracers were not modified by processes that occur before watercloud between commit 633 and this one.

File size: 16.1 KB
Line 
1       SUBROUTINE watercloud(ngrid,nlay,ptimestep,
2     &                pplev,pplay,pdpsrf,pzlay,pt,pdt,
3     &                pq,pdq,pdqcloud,pdtcloud,
4     &                nq,tau,tauscaling,rdust,rice,nuice,
5     &                rsedcloud,rhocloud)
6! to use  'getin'
7      USE ioipsl_getincom
8      IMPLICIT NONE
9
10
11c=======================================================================
12c  Water-ice cloud formation
13
14c  Includes two different schemes:
15c    - A simplified scheme (see simpleclouds.F)
16c    - An improved microphysical scheme (see improvedclouds.F)
17c
18c  There is a time loop specific to cloud formation
19c  due to timescales smaller than the GCM integration timestep.
20c
21c  Authors: Franck Montmessin, Francois Forget, Ehouarn Millour,
22c           J.-B. Madeleine, Thomas Navarro
23c
24c  2004 - 2012
25c=======================================================================
26
27c-----------------------------------------------------------------------
28c   declarations:
29c   -------------
30
31#include "dimensions.h"
32#include "dimphys.h"
33#include "comcstfi.h"
34#include "callkeys.h"
35#include "tracer.h"
36#include "comgeomfi.h"
37#include "dimradmars.h"
38
39c   Inputs:
40c   ------
41
42      INTEGER ngrid,nlay
43      INTEGER nq                 ! nombre de traceurs
44      REAL ptimestep             ! pas de temps physique (s)
45      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
46      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
47      REAL pdpsrf(ngrid)         ! tendence surf pressure
48      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
49      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
50      REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
51
52      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
53      real pdq(ngrid,nlay,nq)    ! tendence avant condensation  (kg/kg.s-1)
54
55      REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
56      REAL tauscaling(ngridmx)   ! Convertion factor for dust amount
57      real rdust(ngridmx,nlay)   ! Dust geometric mean radius (m)
58
59c   Outputs:
60c   -------
61
62      real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
63      REAL pdtcloud(ngrid,nlay)    ! tendence temperature due
64                                   ! a la chaleur latente
65
66      REAL rice(ngrid,nlay)    ! Ice mass mean radius (m)
67                               ! (r_c in montmessin_2004)
68      REAL nuice(ngrid,nlay)   ! Estimated effective variance
69                               !   of the size distribution
70      real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius
71      real rhocloud(ngridmx,nlay)  ! Cloud density (kg.m-3)
72
73c   local:
74c   ------
75     
76      ! for ice radius computation
77      REAL Mo,No
78      REAl ccntyp
79     
80      ! for time loop
81      INTEGER microstep  ! time subsampling step variable
82      INTEGER imicro     ! time subsampling for coupled water microphysics & sedimentation
83      SAVE imicro
84      REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
85      SAVE microtimestep
86     
87      ! tendency given by clouds (inside the micro loop)
88      REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud
89      REAL subpdtcloud(ngrid,nlay)    ! cf. pdtcloud
90
91      ! global tendency (clouds+physics)
92      REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
93      REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
94     
95      REAL CBRT
96      EXTERNAL CBRT
97
98
99      INTEGER iq,ig,l
100      LOGICAL,SAVE :: firstcall=.true.
101
102c    ** un petit test de coherence
103c       --------------------------
104
105      IF (firstcall) THEN
106        IF(ngrid.NE.ngridmx) THEN
107            PRINT*,'STOP dans watercloud'
108            PRINT*,'probleme de dimensions :'
109            PRINT*,'ngrid  =',ngrid
110            PRINT*,'ngridmx  =',ngridmx
111            STOP
112        ENDIF
113         
114        if (nq.gt.nqmx) then
115           write(*,*) 'stop in watercloud (nq.gt.nqmx)!'
116           write(*,*) 'nq=',nq,' nqmx=',nqmx
117           stop
118        endif
119         
120        write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
121        write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice
122               
123        write(*,*) "time subsampling for microphysic ?"
124#ifdef MESOSCALE
125        imicro = 2
126#else
127        imicro = 15
128#endif
129        call getin("imicro",imicro)
130        write(*,*)"imicro = ",imicro
131       
132        microtimestep = ptimestep/real(imicro)
133        write(*,*)"Physical timestep is",ptimestep
134        write(*,*)"Microphysics timestep is",microtimestep
135
136        firstcall=.false.
137      ENDIF ! of IF (firstcall)
138     
139c-----Initialization
140      subpdq(1:ngrid,1:nlay,1:nq) = 0
141      subpdt(1:ngrid,1:nlay)      = 0
142     
143      ! default value if no ice
144      rhocloud(1:ngrid,1:nlay) = rho_dust
145
146
147
148c------------------------------------------------------------------
149c Time subsampling for microphysics
150c------------------------------------------------------------------
151      DO microstep=1,imicro
152     
153c-------------------------------------------------------------------
154c   1.  Tendencies:
155c------------------
156
157
158c------ Temperature tendency subpdt
159        ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
160        ! If imicro=1 subpdt is the same as pdt
161        DO l=1,nlay
162          DO ig=1,ngrid
163             subpdt(ig,l) = subpdt(ig,l)
164     &        + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
165          ENDDO
166        ENDDO
167c------ Tracers tendencies subpdq
168c------ At each micro timestep we add pdq in order to have a stepped entry
169        IF (microphys) THEN
170          DO l=1,nlay
171            DO ig=1,ngrid
172              subpdq(ig,l,igcm_dust_mass) =
173     &            subpdq(ig,l,igcm_dust_mass)
174     &          + pdq(ig,l,igcm_dust_mass)
175              subpdq(ig,l,igcm_dust_number) =
176     &            subpdq(ig,l,igcm_dust_number)
177     &          + pdq(ig,l,igcm_dust_number)
178              subpdq(ig,l,igcm_ccn_mass) =
179     &            subpdq(ig,l,igcm_ccn_mass)
180     &          + pdq(ig,l,igcm_ccn_mass)
181              subpdq(ig,l,igcm_ccn_number) =
182     &            subpdq(ig,l,igcm_ccn_number)
183     &          + pdq(ig,l,igcm_ccn_number)
184            ENDDO
185          ENDDO
186        ENDIF
187        DO l=1,nlay
188          DO ig=1,ngrid
189            subpdq(ig,l,igcm_h2o_ice) =
190     &          subpdq(ig,l,igcm_h2o_ice)
191     &        + pdq(ig,l,igcm_h2o_ice)
192            subpdq(ig,l,igcm_h2o_vap) =
193     &          subpdq(ig,l,igcm_h2o_vap)
194     &        + pdq(ig,l,igcm_h2o_vap)
195          ENDDO
196        ENDDO
197       
198       
199c-------------------------------------------------------------------
200c   2.  Main call to the different cloud schemes:
201c------------------------------------------------
202        IF (microphys) THEN
203           CALL improvedclouds(ngrid,nlay,microtimestep,
204     &             pplay,pt,subpdt,
205     &             pq,subpdq,subpdqcloud,subpdtcloud,
206     &             nq,tauscaling)
207
208        ELSE
209           CALL simpleclouds(ngrid,nlay,microtimestep,
210     &             pplay,pzlay,pt,subpdt,
211     &             pq,subpdq,subpdqcloud,subpdtcloud,
212     &             nq,tau,rice)
213        ENDIF
214
215
216c-------------------------------------------------------------------
217c   3.  Updating tendencies after cloud scheme:
218c-----------------------------------------------
219
220        IF (microphys) THEN
221          DO l=1,nlay
222            DO ig=1,ngrid
223              subpdq(ig,l,igcm_dust_mass) =
224     &            subpdq(ig,l,igcm_dust_mass)
225     &          + subpdqcloud(ig,l,igcm_dust_mass)
226              subpdq(ig,l,igcm_dust_number) =
227     &            subpdq(ig,l,igcm_dust_number)
228     &          + subpdqcloud(ig,l,igcm_dust_number)
229              subpdq(ig,l,igcm_ccn_mass) =
230     &            subpdq(ig,l,igcm_ccn_mass)
231     &          + subpdqcloud(ig,l,igcm_ccn_mass)
232              subpdq(ig,l,igcm_ccn_number) =
233     &            subpdq(ig,l,igcm_ccn_number)
234     &          + subpdqcloud(ig,l,igcm_ccn_number)
235            ENDDO
236          ENDDO
237        ENDIF
238        DO l=1,nlay
239          DO ig=1,ngrid
240            subpdq(ig,l,igcm_h2o_ice) =
241     &          subpdq(ig,l,igcm_h2o_ice)
242     &        + subpdqcloud(ig,l,igcm_h2o_ice)
243            subpdq(ig,l,igcm_h2o_vap) =
244     &          subpdq(ig,l,igcm_h2o_vap)
245     &        + subpdqcloud(ig,l,igcm_h2o_vap)
246          ENDDO
247        ENDDO
248     
249 
250      ENDDO ! of DO microstep=1,imicro
251     
252c-------------------------------------------------------------------
253c   6.  Compute final tendencies after time loop:
254c------------------------------------------------
255
256c------ Temperature tendency pdtcloud
257       DO l=1,nlay
258         DO ig=1,ngrid
259             pdtcloud(ig,l) =
260     &         subpdt(ig,l)/real(imicro)-pdt(ig,l)
261          ENDDO
262       ENDDO
263c------ Tracers tendencies pdqcloud
264       DO l=1,nlay
265         DO ig=1,ngrid
266            pdqcloud(ig,l,igcm_h2o_ice) =
267     &        subpdq(ig,l,igcm_h2o_ice)/real(imicro)
268     &       - pdq(ig,l,igcm_h2o_ice)
269            pdqcloud(ig,l,igcm_h2o_vap) =
270     &        subpdq(ig,l,igcm_h2o_vap)/real(imicro)
271     &       - pdq(ig,l,igcm_h2o_vap)
272          IF(microphys) THEN
273            pdqcloud(ig,l,igcm_dust_mass) =
274     &        subpdq(ig,l,igcm_dust_mass)/real(imicro)
275     &       - pdq(ig,l,igcm_dust_mass)
276            pdqcloud(ig,l,igcm_dust_number) =
277     &        subpdq(ig,l,igcm_dust_number)/real(imicro)
278     &       - pdq(ig,l,igcm_dust_number)
279            pdqcloud(ig,l,igcm_ccn_mass) =
280     &        subpdq(ig,l,igcm_ccn_mass)/real(imicro)
281     &       - pdq(ig,l,igcm_ccn_mass)
282            pdqcloud(ig,l,igcm_ccn_number) =
283     &        subpdq(ig,l,igcm_ccn_number)/real(imicro)
284     &       - pdq(ig,l,igcm_ccn_number)
285           ENDIF
286         ENDDO
287       ENDDO
288
289c------- Due to stepped entry, other processes tendencies can add up to negative values
290c------- Therefore, enforce positive values and conserve mass
291
292       IF(microphys) THEN
293        DO l=1,nlay
294         DO ig=1,ngrid
295          IF ((pq(ig,l,igcm_dust_number) +
296     &      ptimestep* (pdq(ig,l,igcm_dust_number) +
297     &        pdqcloud(ig,l,igcm_dust_number)) .le. 1.)
298     &   .or. (pq(ig,l,igcm_dust_mass) +
299     &      ptimestep* (pdq(ig,l,igcm_dust_mass) +
300     &        pdqcloud(ig,l,igcm_dust_mass)) .le. 1.e-20)) THEN
301         pdqcloud(ig,l,igcm_dust_number) =
302     &     - pq(ig,l,igcm_dust_number)/ptimestep
303     &     - pdq(ig,l,igcm_dust_number) + 1.
304         pdqcloud(ig,l,igcm_ccn_number) = 
305     &     -pdqcloud(ig,l,igcm_dust_number)
306         pdqcloud(ig,l,igcm_dust_mass) =
307     &     - pq(ig,l,igcm_dust_mass)/ptimestep
308     &     - pdq(ig,l,igcm_dust_mass) + 1.e-20
309         pdqcloud(ig,l,igcm_ccn_mass) =
310     &     -pdqcloud(ig,l,igcm_dust_mass)
311          ENDIF
312          IF ((pq(ig,l,igcm_ccn_number) +
313     &      ptimestep* (pdq(ig,l,igcm_ccn_number) +
314     &        pdqcloud(ig,l,igcm_ccn_number)) .le. 1.)
315     &   .or. (pq(ig,l,igcm_ccn_mass) +
316     &      ptimestep* (pdq(ig,l,igcm_ccn_mass) +
317     &        pdqcloud(ig,l,igcm_ccn_mass)) .le. 1.e-20)) THEN
318         pdqcloud(ig,l,igcm_ccn_number) =
319     &     - pq(ig,l,igcm_ccn_number)/ptimestep
320     &     - pdq(ig,l,igcm_ccn_number) + 1.
321         pdqcloud(ig,l,igcm_dust_number) = 
322     &     -pdqcloud(ig,l,igcm_ccn_number)
323         pdqcloud(ig,l,igcm_ccn_mass) =
324     &     - pq(ig,l,igcm_ccn_mass)/ptimestep
325     &     - pdq(ig,l,igcm_ccn_mass) + 1.e-20
326         pdqcloud(ig,l,igcm_dust_mass) =
327     &     -pdqcloud(ig,l,igcm_ccn_mass)
328          ENDIF
329         ENDDO
330        ENDDO
331       ENDIF
332
333        DO l=1,nlay
334         DO ig=1,ngrid
335          IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
336     &       (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice))
337     &       .le. 1.e-8) THEN
338           pdqcloud(ig,l,igcm_h2o_ice) =
339     &     - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
340           pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
341          ENDIF
342          IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
343     &       (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
344     &       .le. 1.e-8) THEN
345           pdqcloud(ig,l,igcm_h2o_vap) =
346     &     - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
347           pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
348          ENDIF
349         ENDDO
350        ENDDO
351
352
353c------Update the ice and dust particle size "rice" for output or photochemistry
354c------Only rsedcloud is used for the water cycle
355
356      IF(scavenging) THEN
357        DO l=1, nlay
358         DO ig=1,ngrid
359
360          rdust(ig,l)=
361     &      CBRT(r3n_q*(pq(ig,l,igcm_dust_mass)+
362     &           (pdq(ig,l,igcm_dust_mass)
363     &           +pdqcloud(ig,l,igcm_dust_mass))*ptimestep)/
364     &      (pq(ig,l,igcm_dust_number)+
365     &           (pdq(ig,l,igcm_dust_number)+
366     &            pdqcloud(ig,l,igcm_dust_number)*ptimestep)))
367          rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
368
369          Mo = pq(ig,l,igcm_h2o_ice)
370     &         + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep
371     &         + (pq(ig,l,igcm_ccn_mass)
372     &          + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
373     &         *tauscaling(ig) + 1.e-30
374          No = (pq(ig,l,igcm_ccn_number)
375     &          + pdqcloud(ig,l,igcm_ccn_number)*ptimestep)
376     &         *tauscaling(ig) + 1.e-30
377          rhocloud(ig,l) = (pq(ig,l,igcm_h2o_ice) + 
378     &                   pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)
379     &                   / Mo * rho_ice
380     &                  + (pq(ig,l,igcm_ccn_mass)
381     &                  + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
382     &                   *tauscaling(ig)/ Mo * rho_dust
383          rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
384          if ((Mo.lt.1.e-15) .or. (No.le.1.)) then
385              rice(ig,l) = 1.e-8
386          else
387              !! AS: only perform computations if conditions not met
388              rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.)
389              rice(ig,l)=min(max(rice(ig,l),1.e-10),500.e-6)
390          endif
391          rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)*
392     &                         (1.+nuice_sed)*(1.+nuice_sed),
393     &                         rdust(ig,l) )
394          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
395         ENDDO
396        ENDDO
397       ELSE
398       
399        if ((igcm_dust_mass.ne.0).and.(igcm_dust_number.ne.0)) then
400         ! recompute rdust(), if we have dust_mass & dust_number tracers
401         DO l=1,nlay
402          DO ig=1,ngrid
403           rdust(ig,l)=
404     &      CBRT(r3n_q*(pq(ig,l,igcm_dust_mass)+
405     &           (pdq(ig,l,igcm_dust_mass)
406     &           +pdqcloud(ig,l,igcm_dust_mass))*ptimestep)/
407     &      (pq(ig,l,igcm_dust_number)+
408     &           (pdq(ig,l,igcm_dust_number)+
409     &            pdqcloud(ig,l,igcm_dust_number)*ptimestep)))
410           rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
411          ENDDO
412         ENDDO
413        endif ! of if ((igcm_dust_mass.ne.0).and.(igcm_dust_number.ne.0))
414       
415        DO l=1,nlay
416          DO ig=1,ngrid
417            ccntyp =
418     &     1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-pzlay(ig,l)/10000.)
419            ccntyp = ccntyp /ccn_factor
420            rice(ig,l)=max( CBRT ( ((pq(ig,l,igcm_h2o_ice)
421     &       + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)/rho_ice
422     &      +ccntyp*(4./3.)*pi*rdust(ig,l)*rdust(ig,l)*rdust(ig,l))
423     &      /(ccntyp*4./3.*pi) ), rdust(ig,l))
424            rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)*
425     &                         (1.+nuice_sed)*(1.+nuice_sed),
426     &                         rdust(ig,l) )
427            rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
428          ENDDO
429        ENDDO
430       ENDIF ! of IF(scavenging)
431       
432       
433! used for rad. transfer calculations
434! nuice is constant because a lognormal distribution is prescribed
435      nuice(1:ngrid,1:nlay)=nuice_ref
436
437
438!--------------------------------------------------------------
439!--------------------------------------------------------------
440     
441
442c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
443c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
444c     Then that should not affect the ice particle radius
445      do ig=1,ngridmx
446        if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
447          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
448     &    rice(ig,2)=rice(ig,3)
449          rice(ig,1)=rice(ig,2)
450        end if
451      end do
452
453c=======================================================================
454
455      END
456 
Note: See TracBrowser for help on using the repository browser.