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

Last change on this file since 1635 was 1467, checked in by tnavarro, 9 years ago

new option supersat to allow for supersaturation of water

File size: 16.3 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      USE comcstfi_h
10      use tracer_mod, only: nqmx, igcm_h2o_vap, igcm_h2o_ice,
11     &                      igcm_dust_mass, igcm_dust_number,
12     &                      igcm_ccn_mass, igcm_ccn_number,
13     &                      rho_dust, nuice_sed, nuice_ref
14      use dimradmars_mod, only: naerkind
15      IMPLICIT NONE
16
17
18c=======================================================================
19c  Water-ice cloud formation
20
21c  Includes two different schemes:
22c    - A simplified scheme (see simpleclouds.F)
23c    - An improved microphysical scheme (see improvedclouds.F)
24c
25c  There is a time loop specific to cloud formation
26c  due to timescales smaller than the GCM integration timestep.
27c
28c  Authors: Franck Montmessin, Francois Forget, Ehouarn Millour,
29c           J.-B. Madeleine, Thomas Navarro
30c
31c  2004 - 2012
32c=======================================================================
33
34c-----------------------------------------------------------------------
35c   declarations:
36c   -------------
37
38#include "callkeys.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(ngrid,naerkind) ! Column dust optical depth at each point
57      REAL tauscaling(ngrid)   ! Convertion factor for dust amount
58      real rdust(ngrid,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(ngrid,nlay) ! Cloud sedimentation radius
72      real rhocloud(ngrid,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      ! no supersaturation when option supersat is false
97      REAL zt(ngrid,nlay)       ! local value of temperature
98      REAL zqsat(ngrid,nlay)    ! saturation
99
100      INTEGER iq,ig,l
101      LOGICAL,SAVE :: firstcall=.true.
102
103c    ** un petit test de coherence
104c       --------------------------
105
106      IF (firstcall) THEN
107         
108        if (nq.gt.nqmx) then
109           write(*,*) 'stop in watercloud (nq.gt.nqmx)!'
110           write(*,*) 'nq=',nq,' nqmx=',nqmx
111           stop
112        endif
113         
114        write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
115        write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice
116               
117        write(*,*) "time subsampling for microphysic ?"
118#ifdef MESOSCALE
119        imicro = 2
120#else
121        imicro = 30
122#endif
123        call getin("imicro",imicro)
124        write(*,*)"imicro = ",imicro
125       
126        microtimestep = ptimestep/real(imicro)
127        write(*,*)"Physical timestep is",ptimestep
128        write(*,*)"Microphysics timestep is",microtimestep
129
130        firstcall=.false.
131      ENDIF ! of IF (firstcall)
132     
133c-----Initialization
134      subpdq(1:ngrid,1:nlay,1:nq) = 0
135      subpdt(1:ngrid,1:nlay)      = 0
136     
137      ! default value if no ice
138      rhocloud(1:ngrid,1:nlay) = rho_dust
139
140
141
142c------------------------------------------------------------------
143c Time subsampling for microphysics
144c------------------------------------------------------------------
145      DO microstep=1,imicro
146     
147c-------------------------------------------------------------------
148c   1.  Tendencies:
149c------------------
150
151
152c------ Temperature tendency subpdt
153        ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
154        ! If imicro=1 subpdt is the same as pdt
155        DO l=1,nlay
156          DO ig=1,ngrid
157             subpdt(ig,l) = subpdt(ig,l)
158     &        + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
159          ENDDO
160        ENDDO
161c------ Tracers tendencies subpdq
162c------ At each micro timestep we add pdq in order to have a stepped entry
163        IF (microphys) THEN
164          DO l=1,nlay
165            DO ig=1,ngrid
166              subpdq(ig,l,igcm_dust_mass) =
167     &            subpdq(ig,l,igcm_dust_mass)
168     &          + pdq(ig,l,igcm_dust_mass)
169              subpdq(ig,l,igcm_dust_number) =
170     &            subpdq(ig,l,igcm_dust_number)
171     &          + pdq(ig,l,igcm_dust_number)
172              subpdq(ig,l,igcm_ccn_mass) =
173     &            subpdq(ig,l,igcm_ccn_mass)
174     &          + pdq(ig,l,igcm_ccn_mass)
175              subpdq(ig,l,igcm_ccn_number) =
176     &            subpdq(ig,l,igcm_ccn_number)
177     &          + pdq(ig,l,igcm_ccn_number)
178            ENDDO
179          ENDDO
180        ENDIF
181        DO l=1,nlay
182          DO ig=1,ngrid
183            subpdq(ig,l,igcm_h2o_ice) =
184     &          subpdq(ig,l,igcm_h2o_ice)
185     &        + pdq(ig,l,igcm_h2o_ice)
186            subpdq(ig,l,igcm_h2o_vap) =
187     &          subpdq(ig,l,igcm_h2o_vap)
188     &        + pdq(ig,l,igcm_h2o_vap)
189          ENDDO
190        ENDDO
191       
192       
193c-------------------------------------------------------------------
194c   2.  Main call to the different cloud schemes:
195c------------------------------------------------
196        IF (microphys) THEN
197           CALL improvedclouds(ngrid,nlay,microtimestep,
198     &             pplay,pt,subpdt,
199     &             pq,subpdq,subpdqcloud,subpdtcloud,
200     &             nq,tauscaling)
201
202        ELSE
203           CALL simpleclouds(ngrid,nlay,microtimestep,
204     &             pplay,pzlay,pt,subpdt,
205     &             pq,subpdq,subpdqcloud,subpdtcloud,
206     &             nq,tau,rice)
207        ENDIF
208
209
210c-------------------------------------------------------------------
211c   3.  Updating tendencies after cloud scheme:
212c-----------------------------------------------
213
214        IF (microphys) THEN
215          DO l=1,nlay
216            DO ig=1,ngrid
217              subpdq(ig,l,igcm_dust_mass) =
218     &            subpdq(ig,l,igcm_dust_mass)
219     &          + subpdqcloud(ig,l,igcm_dust_mass)
220              subpdq(ig,l,igcm_dust_number) =
221     &            subpdq(ig,l,igcm_dust_number)
222     &          + subpdqcloud(ig,l,igcm_dust_number)
223              subpdq(ig,l,igcm_ccn_mass) =
224     &            subpdq(ig,l,igcm_ccn_mass)
225     &          + subpdqcloud(ig,l,igcm_ccn_mass)
226              subpdq(ig,l,igcm_ccn_number) =
227     &            subpdq(ig,l,igcm_ccn_number)
228     &          + subpdqcloud(ig,l,igcm_ccn_number)
229            ENDDO
230          ENDDO
231        ENDIF
232        DO l=1,nlay
233          DO ig=1,ngrid
234            subpdq(ig,l,igcm_h2o_ice) =
235     &          subpdq(ig,l,igcm_h2o_ice)
236     &        + subpdqcloud(ig,l,igcm_h2o_ice)
237            subpdq(ig,l,igcm_h2o_vap) =
238     &          subpdq(ig,l,igcm_h2o_vap)
239     &        + subpdqcloud(ig,l,igcm_h2o_vap)
240          ENDDO
241        ENDDO
242       
243        IF (activice) THEN
244          DO l=1,nlay
245            DO ig=1,ngrid
246              subpdt(ig,l) =
247     &            subpdt(ig,l) + subpdtcloud(ig,l)
248            ENDDO
249          ENDDO
250        ENDIF
251     
252 
253      ENDDO ! of DO microstep=1,imicro
254     
255c-------------------------------------------------------------------
256c   6.  Compute final tendencies after time loop:
257c------------------------------------------------
258
259c------ Temperature tendency pdtcloud
260       DO l=1,nlay
261         DO ig=1,ngrid
262             pdtcloud(ig,l) =
263     &         subpdt(ig,l)/real(imicro)-pdt(ig,l)
264          ENDDO
265       ENDDO
266       
267c------ Tracers tendencies pdqcloud
268       DO l=1,nlay
269         DO ig=1,ngrid
270            pdqcloud(ig,l,igcm_h2o_ice) =
271     &        subpdq(ig,l,igcm_h2o_ice)/real(imicro)
272     &       - pdq(ig,l,igcm_h2o_ice)
273            pdqcloud(ig,l,igcm_h2o_vap) =
274     &        subpdq(ig,l,igcm_h2o_vap)/real(imicro)
275     &       - pdq(ig,l,igcm_h2o_vap)
276         ENDDO
277       ENDDO
278       
279       IF(microphys) THEN
280        DO l=1,nlay
281         DO ig=1,ngrid
282            pdqcloud(ig,l,igcm_ccn_mass) =
283     &        subpdq(ig,l,igcm_ccn_mass)/real(imicro)
284     &       - pdq(ig,l,igcm_ccn_mass)
285            pdqcloud(ig,l,igcm_ccn_number) =
286     &        subpdq(ig,l,igcm_ccn_number)/real(imicro)
287     &       - pdq(ig,l,igcm_ccn_number)
288         ENDDO
289        ENDDO
290       ENDIF
291       
292       IF(scavenging) THEN
293        DO l=1,nlay
294         DO ig=1,ngrid
295            pdqcloud(ig,l,igcm_dust_mass) =
296     &        subpdq(ig,l,igcm_dust_mass)/real(imicro)
297     &       - pdq(ig,l,igcm_dust_mass)
298            pdqcloud(ig,l,igcm_dust_number) =
299     &        subpdq(ig,l,igcm_dust_number)/real(imicro)
300     &       - pdq(ig,l,igcm_dust_number)
301         ENDDO
302        ENDDO
303       ENDIF
304
305c------- Due to stepped entry, other processes tendencies can add up to negative values
306c------- Therefore, enforce positive values and conserve mass
307
308
309       IF(microphys) THEN
310        DO l=1,nlay
311         DO ig=1,ngrid
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       IF(scavenging) THEN
334        DO l=1,nlay
335         DO ig=1,ngrid
336          IF ((pq(ig,l,igcm_dust_number) +
337     &      ptimestep* (pdq(ig,l,igcm_dust_number) +
338     &        pdqcloud(ig,l,igcm_dust_number)) .le. 1.)
339     &   .or. (pq(ig,l,igcm_dust_mass) +
340     &      ptimestep* (pdq(ig,l,igcm_dust_mass) +
341     &        pdqcloud(ig,l,igcm_dust_mass)) .le. 1.e-20)) THEN
342         pdqcloud(ig,l,igcm_dust_number) =
343     &     - pq(ig,l,igcm_dust_number)/ptimestep
344     &     - pdq(ig,l,igcm_dust_number) + 1.
345         pdqcloud(ig,l,igcm_ccn_number) = 
346     &     -pdqcloud(ig,l,igcm_dust_number)
347         pdqcloud(ig,l,igcm_dust_mass) =
348     &     - pq(ig,l,igcm_dust_mass)/ptimestep
349     &     - pdq(ig,l,igcm_dust_mass) + 1.e-20
350         pdqcloud(ig,l,igcm_ccn_mass) =
351     &     -pdqcloud(ig,l,igcm_dust_mass)
352          ENDIF
353         ENDDO
354        ENDDO
355       ENDIF
356
357        DO l=1,nlay
358         DO ig=1,ngrid
359          IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
360     &       (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice))
361     &       .le. 1.e-8) THEN
362           pdqcloud(ig,l,igcm_h2o_ice) =
363     &     - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
364           pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
365          ENDIF
366          IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
367     &       (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
368     &       .le. 1.e-8) THEN
369           pdqcloud(ig,l,igcm_h2o_vap) =
370     &     - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
371           pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
372          ENDIF
373         ENDDO
374        ENDDO
375
376
377c------Update the ice and dust particle size "rice" for output or photochemistry
378c------Only rsedcloud is used for the water cycle
379
380      IF(scavenging) THEN
381        DO l=1, nlay
382         DO ig=1,ngrid
383
384        call updaterdust(
385     &    pq(ig,l,igcm_dust_mass) +                   ! dust mass
386     &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
387     &    pdqcloud(ig,l,igcm_dust_mass))*ptimestep,   ! dust mass
388     &    pq(ig,l,igcm_dust_number) +                 ! dust number
389     &   (pdq(ig,l,igcm_dust_number) +                ! dust number
390     &    pdqcloud(ig,l,igcm_dust_number))*ptimestep, ! dust number
391     &    rdust(ig,l))
392
393         ENDDO
394        ENDDO
395      ENDIF
396
397      IF(microphys) THEN
398
399       ! In case one does not want to allow supersatured water when using microphysics.
400       ! Not done by default.
401       IF(.not.supersat) THEN     
402        zt  = pt + (pdt+pdtcloud)*ptimestep
403        call watersat(ngrid*nlay,zt,pplay,zqsat)
404        DO l=1, nlay
405         DO ig=1,ngrid
406          IF (pq(ig,l,igcm_h2o_vap)
407     &      + (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
408     &      * ptimestep .ge. zqsat(ig,l)) THEN
409             pdqcloud(ig,l,igcm_h2o_vap) =
410     &         (zqsat(ig,l) - pq(ig,l,igcm_h2o_vap))/ptimestep
411     &        - pdq(ig,l,igcm_h2o_vap)
412             pdqcloud(ig,l,igcm_h2o_ice) =
413     &         -pdqcloud(ig,l,igcm_h2o_vap)
414             ! no need to correct ccn_number, updaterad can handle this properly.
415          ENDIF
416         ENDDO
417        ENDDO       
418       ENDIF
419       
420       DO l=1, nlay
421         DO ig=1,ngrid
422
423        call updaterice_micro(
424     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
425     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
426     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
427     &    pq(ig,l,igcm_ccn_mass) +                   ! ccn mass
428     &   (pdq(ig,l,igcm_ccn_mass) +                  ! ccn mass
429     &    pdqcloud(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
430     &    pq(ig,l,igcm_ccn_number) +                 ! ccn number
431     &   (pdq(ig,l,igcm_ccn_number) +                ! ccn number
432     &    pdqcloud(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
433     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
434         
435         ENDDO
436       ENDDO
437       
438      ELSE ! no microphys
439       
440        DO l=1,nlay
441          DO ig=1,ngrid
442         
443        call updaterice_typ(
444     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
445     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
446     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
447     &    tau(ig,1),pzlay(ig,l),rice(ig,l))
448
449          ENDDO
450         ENDDO
451       
452       ENDIF ! of IF(microphys)
453     
454     
455     
456c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
457c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
458c     Then that should not affect the ice particle radius
459      do ig=1,ngrid
460        if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
461          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
462     &    rice(ig,2)=rice(ig,3)
463          rice(ig,1)=rice(ig,2)
464        end if
465      end do
466       
467       
468       DO l=1,nlay
469         DO ig=1,ngrid
470           rsedcloud(ig,l)=max(rice(ig,l)*
471     &                 (1.+nuice_sed)*(1.+nuice_sed)*(1.+nuice_sed),
472     &                    rdust(ig,l))
473!          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
474         ENDDO
475       ENDDO
476       
477! used for rad. transfer calculations
478! nuice is constant because a lognormal distribution is prescribed
479      nuice(1:ngrid,1:nlay)=nuice_ref
480
481
482
483c=======================================================================
484
485      END
486 
Note: See TracBrowser for help on using the repository browser.