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

Last change on this file since 937 was 882, checked in by tnavarro, 12 years ago

small correction in makegcm, latent heat in clouds, changed threshold in dust radius

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 = 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        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.