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

Last change on this file since 823 was 746, checked in by jbmadeleine, 12 years ago

Changed tau(ig,l) to tau(ig,1) in the call to updaterice_typ in
watercloud.F (resulted in NaNs?, we were out of the bounds of tau).

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