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

Last change on this file since 342 was 334, checked in by aslmd, 14 years ago

LMDZ.MARS:

28/10/11 == FL + AS

ADDED THE NEW FRAMEWORK FOR PHOTOCHEMISTRY
This is not the last version. A new rewritten version of calchim.F [using LAPACK] is yet to be committed.
--> A new version for photochemistry routines : *_B.F no longer exists, new routines in aeronomars
D 333 libf/aeronomars/photochemist_B.F
D 333 libf/aeronomars/init_chimie_B.F
A 0 libf/aeronomars/read_phototable.F
M 333 libf/aeronomars/calchim.F
A 0 libf/aeronomars/photochemistry.F
M 333 libf/aeronomars/chimiedata.h
--> Changed calls to calchim and watercloud [surfdust and surfice needed] in physiq.F
--> Left commented code for outputs in physiq.F [search for FL]
--> Added settings which works for 35 levels in inidissip.F according to FL. Commented for the moment.
--> Checked compilation and run, looks fine. Note that 'jmars.20101220' is needed.

File size: 13.0 KB
Line 
1       SUBROUTINE watercloud(ngrid,nlay, ptimestep,
2     &                pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
3     &                pq,pdq,pdqcloud,pdqscloud,pdtcloud,
4     &                nq,naersize,tau,
5     &                ccn,rdust,rice,nuice,
6     &                surfdust,surfice)
7      IMPLICIT NONE
8
9c=======================================================================
10c     Treatment of saturation of water vapor
11c
12c
13c     Modif de zq si saturation dans l'atmosphere
14c     si zq(ig,l)> zqsat(ig,l) ->    zq(ig,l)=zqsat(ig,l)
15c     Le test est effectue de bas en haut. L'eau condensee
16c    (si saturation) est remise dans la couche en dessous.
17c     L'eau condensee dans la couche du bas est deposee a la surface
18c       
19c    Modification: Franck Montmessin water ice scheme
20c                  Francois Forget : change nuclei density & outputs   
21c                  Ehouarn Millour: sept.2008, tracers are now handled
22c                                   by name (and not fixed index)
23c
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
37c   Inputs:
38c   ------
39
40      INTEGER ngrid,nlay
41      REAL ptimestep             ! pas de temps physique (s)
42      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
43      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
44      REAL pdpsrf(ngrid)         ! tendance surf pressure
45      REAL pzlev(ngrid,nlay+1)   ! altitude at layer boundaries
46      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
47      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
48      REAL pdt(ngrid,nlay)       ! tendance temperature des autres param.
49
50      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
51      real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
52
53      integer nq         ! nombre de traceurs
54      integer naersize   ! nombre de traceurs radiativement actifs (=naerkind)
55      REAL tau(ngridmx,naersize)
56      REAL ccn(ngridmx,nlayermx)   ! Cloud condensation nuclei
57                                   !   (particules kg-1)
58      real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
59
60c   Outputs:
61c   -------
62
63      real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation H2O(kg/kg.s-1)
64      real pdqscloud(ngrid,nq)     ! flux en surface (kg.m-2.s-1)
65      REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
66                                   !   a la chaleur latente
67
68      REAL rice(ngrid,nlay)    ! Ice mass mean radius (m)
69                               ! (r_c in montmessin_2004)
70      REAL nuice(ngrid,nlay)   ! Estimated effective variance
71                               !   of the size distribution
72
73      REAL surfdust(ngrid,nlay)   ! dust surface area (micron2/cm3, if photochemistry)
74      REAL surfice(ngrid,nlay)    !  ice surface area (micron2/cm3, if photochemistry)
75
76c   local:
77c   ------
78
79      REAL CBRT
80      EXTERNAL CBRT
81      INTEGER ig,l
82
83
84      REAL zq(ngridmx,nlayermx,nqmx)  ! local value of tracers
85      REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
86      REAL zqsat(ngridmx,nlayermx)    ! saturation
87      REAL zt(ngridmx,nlayermx)       ! local value of temperature
88
89      REAL masse (ngridmx,nlayermx)
90      REAL epaisseur (ngridmx,nlayermx)
91      REAL npart(ngridmx,nlayermx)    !  Cloud condensation nuclei (#.cm-3)
92!      REAL rfinal        ! Ice crystal radius after condensation(m)
93      REAL*8 seq           ! Equilibrium saturation ration (accounting for curvature effect)
94      REAL*8 dzq           ! masse de glace echangee (kg/kg)
95      REAL lw       !Latent heat of sublimation (J.kg-1)
96      REAL,PARAMETER :: To=273.15 ! reference temperature, T=273.15 K
97
98      REAL Ctot
99      REAL*8 ph2o,satu
100      REAL*8 gr,Cste,up,dwn,newvap
101
102      LOGICAL,SAVE :: firstcall=.true.
103! To use more refined microphysics, set improved to .true.
104      LOGICAL,PARAMETER :: improved=.true.
105
106c     Pour diagnostique :
107c     ~~~~~~~~~~~~~~~~~
108c     REAL icetot(ngridmx)             ! Total mass of water ice (kg/m2)
109c     REAL rave(ngridmx)               ! Mean crystal radius in a column (m)
110
111! indexes of water vapour, water ice and dust tracers:
112      INTEGER,SAVE :: i_h2o=0  ! water vapour
113      INTEGER,SAVE :: i_ice=0  ! water ice
114
115c    ** un petit test de coherence
116c       --------------------------
117
118      IF (firstcall) THEN
119        IF(ngrid.NE.ngridmx) THEN
120            PRINT*,'STOP dans watercloud'
121            PRINT*,'probleme de dimensions :'
122            PRINT*,'ngrid  =',ngrid
123            PRINT*,'ngridmx  =',ngridmx
124            STOP
125        ENDIF
126         
127        if (nq.gt.nqmx) then
128           write(*,*) 'stop in watercloud (nq.gt.nqmx)!'
129           write(*,*) 'nq=',nq,' nqmx=',nqmx
130           stop
131        endif
132         
133        i_h2o=igcm_h2o_vap
134        i_ice=igcm_h2o_ice
135       
136        write(*,*) "watercloud: i_h2o=",i_h2o
137        write(*,*) "            i_ice=",i_ice
138
139        firstcall=.false.
140      ENDIF ! of IF (firstcall)
141
142
143c-----------------------------------------------------------------------
144c    1. initialisation
145c    -----------------
146
147c    On "update" la valeur de q(nqmx) (water vapor) et temperature.
148c    On effectue qqes calculs preliminaires sur les couches :
149c    masse (kg.m-2), epaisseur(m).
150
151      do l=1,nlay
152        do ig=1,ngrid
153          zq(ig,l,i_h2o)=pq(ig,l,i_h2o)+pdq(ig,l,i_h2o)*ptimestep
154          zq(ig,l,i_h2o)=max(zq(ig,l,i_h2o),1.E-30) ! FF 12/2004
155          zq0(ig,l,i_h2o)=zq(ig,l,i_h2o)
156          zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
157          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
158          epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
159
160          zq(ig,l,i_ice)=pq(ig,l,i_ice)+pdq(ig,l,i_ice)*ptimestep
161          zq(ig,l,i_ice)=max(zq(ig,l,i_ice),0.) ! FF 12/2004
162          zq0(ig,l,i_ice)=zq(ig,l,i_ice)
163
164c         This typical profile is not used anymore; rdust is
165c           set up in updatereffrad.F.
166c         rdust(ig,l)= max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9)
167        enddo
168      enddo
169
170      do l=1,nlay
171        do ig=1,ngrid
172c         Calcul du rayon moyen des particules de glace.
173c         Hypothese : Dans une couche, la glace presente se
174c         repartit uniformement autour du nbre de poussieres
175c         dont le rayon moyen est prescrit par rdust.
176c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
177          rice(ig,l)=CBRT( ( zq(ig,l,i_ice)/rho_ice+
178     &       ccn(ig,l)*(4./3.)*pi*rdust(ig,l)**3. )
179     &       / (ccn(ig,l)*4./3.*pi) )
180          rice(ig,l)=max(rice(ig,l),rdust(ig,l))
181c         Effective variance of the size distribution
182          nuice(ig,l)=nuice_ref
183        enddo ! of do ig=1,ngrid
184      enddo ! of dol=1,nlay
185
186      if (photochem) then
187c        computation of dust and ice surface area (micron2/cm3)
188c        for heterogeneous chemistry
189
190         do l = 1,nlay
191            do ig = 1,ngrid
192c                                                                     
193c              npart: number density of ccn in #/cm3   
194c                                                     
195               npart(ig,l) = 1.e-6*ccn(ig,l) 
196     $                       *masse(ig,l)/epaisseur(ig,l)
197c
198c              dust and ice surface area
199c                                                                 
200               surfdust(ig,l) = npart(ig,l)*4.*pi*1.e12*rdust(ig,l)**2
201c                                                                 
202               if (rice(ig,l) .ge. rdust(ig,l)) then                   
203                  surfice(ig,l)  = npart(ig,l)*4.*pi*1.e12*rice(ig,l)**2
204                  surfdust(ig,l) = 0.
205               else                                                   
206                  surfice(ig,l) = 0.                                 
207               end if                                             
208            end do      ! of do ig=1,ngrid
209         end do         ! of do l=1,nlay
210      end if            ! of photochem
211c
212      pdqscloud(1:ngrid,1:nq)=0
213      pdqcloud(1:ngrid,1:nlay,1:nq)=0
214      pdtcloud(1:ngrid,1:nlay)=0
215
216c     icetot(1:ngrid)=0
217c     rave(1:ngrid)=0
218
219c    ----------------------------------------------
220c
221c
222c       Rapport de melange a saturation dans la couche l : -------
223c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
224
225        call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
226
227c       taux de condensation (kg/kg/s-1) dans les differentes couches
228c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
229
230c       Iceparty is not used anymore: water=>iceparty (JBM).
231c       if(iceparty) then
232
233        do l=1,nlay
234          do ig=1,ngrid
235
236          IF (improved) then
237c         Improved microphysics scheme
238c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
239           Ctot = zq(ig,l,i_h2o) + zq(ig,l,i_ice)
240           ph2o = zq(ig,l,i_h2o) * 44. / 18. * pplay(ig,l)
241           satu = zq(ig,l,i_h2o) / zqsat(ig,l)
242
243           call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
244     &        ph2o,ph2o/satu,seq,rice(ig,l),gr)
245           Cste = ptimestep * 4. * pi * rice(ig,l)
246     *              * rho_ice * ccn(ig,l)
247           up   = zq(ig,l,i_h2o) + Cste * gr * seq
248           dwn  =    1.       + Cste * gr / zqsat(ig,l)
249           newvap = min(up/dwn,Ctot)
250
251           gr  = gr * ( newvap/zqsat(ig,l) - seq )
252           dzq = min( max( Cste * gr,-zq(ig,l,i_ice) )
253     *             , zq(ig,l,i_h2o) )
254
255c          Nucleation (sat ratio must be larger than a critical value)
256c          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
257           if (satu.gt.1.) then
258             if (satu.le.1.4.and.zq(ig,l,i_ice).lt.1.e-8)
259     *          dzq = 0.
260           endif
261
262           ELSE
263c          Old version
264c          ~~~~~~~~~~~
265           if (zq(ig,l,i_h2o).ge.zqsat(ig,l))then  !  Condensation
266             dzq=zq(ig,l,i_h2o)-zqsat(ig,l)               
267           elseif(zq(ig,l,i_h2o).lt.zqsat(ig,l))then  ! Sublimation
268             dzq=-min(zqsat(ig,l)-zq(ig,l,i_h2o),zq(ig,l,i_ice))
269           endif
270
271           ENDIF ! of IF (improved)
272
273c           Water Mass change
274c           ~~~~~~~~~~~~~~~~~
275            zq(ig,l,i_ice)=zq(ig,l,i_ice)+dzq
276            zq(ig,l,i_h2o)=zq(ig,l,i_h2o)-dzq
277
278            rice(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_ice
279     &       +ccn(ig,l)*(4./3.)*pi*rdust(ig,l)**3.)
280     &       /(ccn(ig,l)*4./3.*pi) ), rdust(ig,l))
281
282            enddo ! of do ig=1,ngrid
283          enddo ! of do l=1,nlay
284
285c       The following part have been commented because iceparty
286c           is not used anymore: water=>iceparty (JBM).
287
288c       else   ! if not iceparty
289
290c         Saturation couche nlay a 2 :
291c         ~~~~~~~~~~~~~~~~~~~~~~~~~~
292c         do l=nlay,2, -1
293c          do ig=1,ngrid
294c           if (zq(ig,l,i_h2o).gt.zqsat(ig,l))then
295c             zq(ig,l-1,i_h2o)= zq(ig,l-1,i_h2o)+
296c    &                          (zq(ig,l,i_h2o)-zqsat(ig,l))
297c    &          *(pplev(ig,l)-pplev(ig,l+1))/(pplev(ig,l-1)-pplev(ig,l))
298c             zq(ig,l,i_h2o)=zqsat(ig,l)
299c           endif
300c          enddo
301c         enddo
302
303c       Saturation couche l=1 si pas iceparty
304c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
305c         do ig=1,ngridmx
306c           if (zq(ig,1,i_h2o).gt.zqsat(ig,1))then
307c             pdqscloud(ig,i_ice)=(zq(ig,1,i_h2o)-zqsat(ig,1))
308c    &           *(pplev(ig,1)-pplev(ig,2))/(g*ptimestep)
309c             zq(ig,1,i_h2o)=zqsat(ig,1)
310c           endif
311c         enddo
312
313c       endif   ! of if (iceparty)
314
315c       Tendance finale
316c       ~~~~~~~~~~~~~~~
317        do l=1, nlay
318          do ig=1,ngridmx
319            pdqcloud(ig,l,i_h2o)=(zq(ig,l,i_h2o)
320     &                              -zq0(ig,l,i_h2o))/ptimestep
321            pdqcloud(ig,l,i_ice) =
322     &        (zq(ig,l,i_ice) - zq0(ig,l,i_ice))/ptimestep
323            lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
324            pdtcloud(ig,l)=-pdqcloud(ig,l,i_h2o)*lw/cpp
325          end do
326        end do
327
328c       A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
329c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
330c       Then that should not affect the ice particle radius
331        do ig=1,ngridmx
332          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
333            if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
334     &      rice(ig,2)=rice(ig,3)
335            rice(ig,1)=rice(ig,2)
336          end if
337        end do
338       
339c**************************************************
340c       Output 
341c**************************************************
342! NB: for diagnostics use zq(), the updated value of tracers
343
344c        do ig=1,ngridmx
345c         do l=1 ,nlay
346c           masse de glace d'eau dans la couche l
347c           icetot(ig)=icetot(ig)+masse(ig,l)*zq(ig,l,i_ice)
348c           rayon moyen des cristaux dans la colonne ig
349c           rave(ig)=rave(ig)+masse(ig,l)*zq(ig,l,i_ice)*rice(ig,l)
350c         enddo
351c         rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
352c         if (icetot(ig)*1000.lt.0.01) rave(ig)=0.
353c        enddo   ! (ngridmx)
354c**************************************************
355
356      RETURN
357      END
358 
Note: See TracBrowser for help on using the repository browser.