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

Last change on this file since 1009 was 951, checked in by tnavarro, 12 years ago

default microphysics timestep changed

File size: 15.4 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 = 30
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        IF (activice) THEN
248          DO l=1,nlay
249            DO ig=1,ngrid
250              subpdt(ig,l) =
251     &            subpdt(ig,l) + subpdtcloud(ig,l)
252            ENDDO
253          ENDDO
254        ENDIF
255     
256 
257      ENDDO ! of DO microstep=1,imicro
258     
259c-------------------------------------------------------------------
260c   6.  Compute final tendencies after time loop:
261c------------------------------------------------
262
263c------ Temperature tendency pdtcloud
264       DO l=1,nlay
265         DO ig=1,ngrid
266             pdtcloud(ig,l) =
267     &         subpdt(ig,l)/real(imicro)-pdt(ig,l)
268          ENDDO
269       ENDDO
270       
271c------ Tracers tendencies pdqcloud
272       DO l=1,nlay
273         DO ig=1,ngrid
274            pdqcloud(ig,l,igcm_h2o_ice) =
275     &        subpdq(ig,l,igcm_h2o_ice)/real(imicro)
276     &       - pdq(ig,l,igcm_h2o_ice)
277            pdqcloud(ig,l,igcm_h2o_vap) =
278     &        subpdq(ig,l,igcm_h2o_vap)/real(imicro)
279     &       - pdq(ig,l,igcm_h2o_vap)
280         ENDDO
281       ENDDO
282       
283       IF(microphys) THEN
284        DO l=1,nlay
285         DO ig=1,ngrid
286            pdqcloud(ig,l,igcm_ccn_mass) =
287     &        subpdq(ig,l,igcm_ccn_mass)/real(imicro)
288     &       - pdq(ig,l,igcm_ccn_mass)
289            pdqcloud(ig,l,igcm_ccn_number) =
290     &        subpdq(ig,l,igcm_ccn_number)/real(imicro)
291     &       - pdq(ig,l,igcm_ccn_number)
292         ENDDO
293        ENDDO
294       ENDIF
295       
296       IF(scavenging) THEN
297        DO l=1,nlay
298         DO ig=1,ngrid
299            pdqcloud(ig,l,igcm_dust_mass) =
300     &        subpdq(ig,l,igcm_dust_mass)/real(imicro)
301     &       - pdq(ig,l,igcm_dust_mass)
302            pdqcloud(ig,l,igcm_dust_number) =
303     &        subpdq(ig,l,igcm_dust_number)/real(imicro)
304     &       - pdq(ig,l,igcm_dust_number)
305         ENDDO
306        ENDDO
307       ENDIF
308
309c------- Due to stepped entry, other processes tendencies can add up to negative values
310c------- Therefore, enforce positive values and conserve mass
311
312
313       IF(microphys) THEN
314        DO l=1,nlay
315         DO ig=1,ngrid
316          IF ((pq(ig,l,igcm_ccn_number) +
317     &      ptimestep* (pdq(ig,l,igcm_ccn_number) +
318     &        pdqcloud(ig,l,igcm_ccn_number)) .le. 1.)
319     &   .or. (pq(ig,l,igcm_ccn_mass) +
320     &      ptimestep* (pdq(ig,l,igcm_ccn_mass) +
321     &        pdqcloud(ig,l,igcm_ccn_mass)) .le. 1.e-20)) THEN
322         pdqcloud(ig,l,igcm_ccn_number) =
323     &     - pq(ig,l,igcm_ccn_number)/ptimestep
324     &     - pdq(ig,l,igcm_ccn_number) + 1.
325         pdqcloud(ig,l,igcm_dust_number) = 
326     &     -pdqcloud(ig,l,igcm_ccn_number)
327         pdqcloud(ig,l,igcm_ccn_mass) =
328     &     - pq(ig,l,igcm_ccn_mass)/ptimestep
329     &     - pdq(ig,l,igcm_ccn_mass) + 1.e-20
330         pdqcloud(ig,l,igcm_dust_mass) =
331     &     -pdqcloud(ig,l,igcm_ccn_mass)
332          ENDIF
333         ENDDO
334        ENDDO
335       ENDIF
336
337       IF(scavenging) THEN
338        DO l=1,nlay
339         DO ig=1,ngrid
340          IF ((pq(ig,l,igcm_dust_number) +
341     &      ptimestep* (pdq(ig,l,igcm_dust_number) +
342     &        pdqcloud(ig,l,igcm_dust_number)) .le. 1.)
343     &   .or. (pq(ig,l,igcm_dust_mass) +
344     &      ptimestep* (pdq(ig,l,igcm_dust_mass) +
345     &        pdqcloud(ig,l,igcm_dust_mass)) .le. 1.e-20)) THEN
346         pdqcloud(ig,l,igcm_dust_number) =
347     &     - pq(ig,l,igcm_dust_number)/ptimestep
348     &     - pdq(ig,l,igcm_dust_number) + 1.
349         pdqcloud(ig,l,igcm_ccn_number) = 
350     &     -pdqcloud(ig,l,igcm_dust_number)
351         pdqcloud(ig,l,igcm_dust_mass) =
352     &     - pq(ig,l,igcm_dust_mass)/ptimestep
353     &     - pdq(ig,l,igcm_dust_mass) + 1.e-20
354         pdqcloud(ig,l,igcm_ccn_mass) =
355     &     -pdqcloud(ig,l,igcm_dust_mass)
356          ENDIF
357         ENDDO
358        ENDDO
359       ENDIF
360
361        DO l=1,nlay
362         DO ig=1,ngrid
363          IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
364     &       (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice))
365     &       .le. 1.e-8) THEN
366           pdqcloud(ig,l,igcm_h2o_ice) =
367     &     - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
368           pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
369          ENDIF
370          IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
371     &       (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
372     &       .le. 1.e-8) THEN
373           pdqcloud(ig,l,igcm_h2o_vap) =
374     &     - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
375           pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
376          ENDIF
377         ENDDO
378        ENDDO
379
380
381c------Update the ice and dust particle size "rice" for output or photochemistry
382c------Only rsedcloud is used for the water cycle
383
384      IF(scavenging) THEN
385        DO l=1, nlay
386         DO ig=1,ngrid
387
388        call updaterdust(
389     &    pq(ig,l,igcm_dust_mass) +                   ! dust mass
390     &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
391     &    pdqcloud(ig,l,igcm_dust_mass))*ptimestep,   ! dust mass
392     &    pq(ig,l,igcm_dust_number) +                 ! dust number
393     &   (pdq(ig,l,igcm_dust_number) +                ! dust number
394     &    pdqcloud(ig,l,igcm_dust_number))*ptimestep, ! dust number
395     &    rdust(ig,l))
396
397         ENDDO
398        ENDDO
399      ENDIF
400       
401       
402      IF(microphys) THEN
403       
404       DO l=1, nlay
405         DO ig=1,ngrid
406
407        call updaterice_micro(
408     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
409     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
410     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
411     &    pq(ig,l,igcm_ccn_mass) +                   ! ccn mass
412     &   (pdq(ig,l,igcm_ccn_mass) +                  ! ccn mass
413     &    pdqcloud(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
414     &    pq(ig,l,igcm_ccn_number) +                 ! ccn number
415     &   (pdq(ig,l,igcm_ccn_number) +                ! ccn number
416     &    pdqcloud(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
417     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
418         
419         ENDDO
420       ENDDO
421       
422      ELSE ! no microphys
423       
424        DO l=1,nlay
425          DO ig=1,ngrid
426         
427        call updaterice_typ(
428     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
429     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
430     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
431     &    tau(ig,1),pzlay(ig,l),rice(ig,l))
432
433          ENDDO
434         ENDDO
435       
436       ENDIF ! of IF(microphys)
437     
438     
439     
440c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
441c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
442c     Then that should not affect the ice particle radius
443      do ig=1,ngridmx
444        if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
445          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
446     &    rice(ig,2)=rice(ig,3)
447          rice(ig,1)=rice(ig,2)
448        end if
449      end do
450       
451       
452       DO l=1,nlay
453         DO ig=1,ngrid
454           rsedcloud(ig,l)=max(rice(ig,l)*
455     &                 (1.+nuice_sed)*(1.+nuice_sed)*(1.+nuice_sed),
456     &                    rdust(ig,l))
457!          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
458         ENDDO
459       ENDDO
460       
461! used for rad. transfer calculations
462! nuice is constant because a lognormal distribution is prescribed
463      nuice(1:ngrid,1:nlay)=nuice_ref
464
465
466
467c=======================================================================
468
469      END
470 
Note: See TracBrowser for help on using the repository browser.