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

Last change on this file since 1242 was 1226, checked in by aslmd, 11 years ago

LMDZ.MARS : Replaced comcstfi and planete includes by modules.

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