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

Last change on this file since 524 was 522, checked in by tnavarro, 13 years ago

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