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

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

default microphysics timestep changed

File size: 15.4 KB
RevLine 
[633]1       SUBROUTINE watercloud(ngrid,nlay,ptimestep,
2     &                pplev,pplay,pdpsrf,pzlay,pt,pdt,
[626]3     &                pq,pdq,pdqcloud,pdtcloud,
[358]4     &                nq,tau,tauscaling,rdust,rice,nuice,
5     &                rsedcloud,rhocloud)
[633]6! to use  'getin'
7      USE ioipsl_getincom
[740]8      USE updaterad
[38]9      IMPLICIT NONE
10
[633]11
[38]12c=======================================================================
[358]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)
[38]18c
[633]19c  There is a time loop specific to cloud formation
20c  due to timescales smaller than the GCM integration timestep.
21c
[358]22c  Authors: Franck Montmessin, Francois Forget, Ehouarn Millour,
[522]23c           J.-B. Madeleine, Thomas Navarro
[38]24c
[633]25c  2004 - 2012
[38]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"
[358]38#include "dimradmars.h"
[38]39
40c   Inputs:
41c   ------
42
43      INTEGER ngrid,nlay
[633]44      INTEGER nq                 ! nombre de traceurs
[38]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)
[633]48      REAL pdpsrf(ngrid)         ! tendence surf pressure
[38]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)
[633]51      REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
[38]52
53      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
[633]54      real pdq(ngrid,nlay,nq)    ! tendence avant condensation  (kg/kg.s-1)
[38]55
[633]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)
[38]59
60c   Outputs:
61c   -------
62
[633]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
[38]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
[633]71      real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius
72      real rhocloud(ngridmx,nlay)  ! Cloud density (kg.m-3)
[38]73
74c   local:
75c   ------
[633]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
[38]91
[633]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
[38]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         
[358]118        write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
119        write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice
[633]120               
121        write(*,*) "time subsampling for microphysic ?"
122#ifdef MESOSCALE
123        imicro = 2
124#else
[951]125        imicro = 30
[633]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
[38]133
134        firstcall=.false.
135      ENDIF ! of IF (firstcall)
[522]136     
[633]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
[38]143
[633]144
145
146c------------------------------------------------------------------
147c Time subsampling for microphysics
148c------------------------------------------------------------------
149      DO microstep=1,imicro
[522]150     
[633]151c-------------------------------------------------------------------
152c   1.  Tendencies:
153c------------------
[38]154
[633]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,
[645]210     &             nq,tau,rice)
[633]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
[882]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
[633]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
[740]270       
[633]271c------ Tracers tendencies pdqcloud
[703]272       DO l=1,nlay
[633]273         DO ig=1,ngrid
[703]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)
[740]280         ENDDO
281       ENDDO
282       
283       IF(microphys) THEN
284        DO l=1,nlay
285         DO ig=1,ngrid
[703]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)
[633]292         ENDDO
[740]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
[633]308
309c------- Due to stepped entry, other processes tendencies can add up to negative values
310c------- Therefore, enforce positive values and conserve mass
311
[740]312
[633]313       IF(microphys) THEN
314        DO l=1,nlay
315         DO ig=1,ngrid
[654]316          IF ((pq(ig,l,igcm_ccn_number) +
[633]317     &      ptimestep* (pdq(ig,l,igcm_ccn_number) +
[654]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
[633]322         pdqcloud(ig,l,igcm_ccn_number) =
323     &     - pq(ig,l,igcm_ccn_number)/ptimestep
[654]324     &     - pdq(ig,l,igcm_ccn_number) + 1.
[633]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
[654]329     &     - pdq(ig,l,igcm_ccn_mass) + 1.e-20
[633]330         pdqcloud(ig,l,igcm_dust_mass) =
331     &     -pdqcloud(ig,l,igcm_ccn_mass)
332          ENDIF
333         ENDDO
334        ENDDO
335       ENDIF
336
[740]337       IF(scavenging) THEN
[633]338        DO l=1,nlay
339         DO ig=1,ngrid
[740]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
[633]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
[740]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))
[633]396
397         ENDDO
398        ENDDO
[740]399      ENDIF
[645]400       
[740]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         
[645]419         ENDDO
[740]420       ENDDO
[645]421       
[740]422      ELSE ! no microphys
423       
[645]424        DO l=1,nlay
425          DO ig=1,ngrid
[740]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
[746]431     &    tau(ig,1),pzlay(ig,l),rice(ig,l))
[740]432
[633]433          ENDDO
[740]434         ENDDO
[633]435       
[740]436       ENDIF ! of IF(microphys)
[633]437     
[740]438     
439     
[358]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
[740]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
[38]464
[740]465
466
[633]467c=======================================================================
468
[38]469      END
470 
Note: See TracBrowser for help on using the repository browser.