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

Last change on this file since 1242 was 57, checked in by aslmd, 14 years ago

mineur LMD_MM_MARS: ajout du GCM ancienne physique, systeme maintenant complet sur SVN (ne manque que la base de donnees d'etats initiaux)

File size: 13.9 KB
Line 
1       SUBROUTINE watercloud(ngrid,nlay, ptimestep,
2     &      pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
3     &                pq,pdq,pdqcloud,qsurf,pdqscloud,pdtcloud,
4     &                nq,nsize,tau,icount,zls)
5      IMPLICIT NONE
6
7c=======================================================================
8c     Treatment of saturation of water vapor
9c
10c
11c     Modif de zq si saturation dans l'atmopshere
12c     si zq(ig,l)> zqsat(ig,l) ->    zq(ig,l)=zqsat(ig,l)
13c     Le test est effectue de bas en haut. L'eau condensee
14c    (si saturation) est remise dans la couche en dessous.
15c     L'eau condensee dans la couche du bas est deposee a la surface
16c
17c     WARNING : water vapor mixing ratio is assumed to be
18c     q(nqmx)
19c       
20c    Modification: Franck Montmessin water ice scheme
21c                  Francois Forget : change nuclei density & outputs   
22c
23c=======================================================================
24
25c-----------------------------------------------------------------------
26c   declarations:
27c   -------------
28
29#include "dimensions.h"
30#include "dimphys.h"
31#include "comcstfi.h"
32#include "callkeys.h"
33#include "tracer.h"
34
35#include "fisice.h"
36#include "comgeomfi.h"
37c
38c   arguments:
39c   ----------
40
41      INTEGER ngrid,nlay,nsize
42      INTEGER icount
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 pzlev(ngrid,nlay+1)   ! altitude at layer boundaries
47      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
48      REAL pdpsrf(ngrid)         ! tendance surf pressure
49      REAL pt(ngrid,nlay)        ! temperature au centre des couches (K)
50      REAL pdt(ngrid,nlay)       ! tendance temperature des autres param.
51      REAL pdtcloud(ngrid,nlay)  ! tendance temperature due a la chaleur latente
52      REAL tau(ngridmx,nsize)
53      REAL zls
54
55c    Traceurs :
56      integer nq                   ! nombre de traceurs
57      real pq(ngrid,nlay,nq)       ! traceur (kg/kg)
58      real pdq(ngrid,nlay,nq)      ! tendance avant condensation  (kg/kg.s-1)
59      real pdqcloud(ngrid,nlay,nq) ! tendance de la condesation H2O(kg/kg.s-1)
60      real qsurf(ngrid,nq)
61      real pdqscloud(ngrid,nq)     ! flux en surface (kg.m-2.s-1)
62     
63c   local:
64c   ------
65
66      REAL CBRT
67      EXTERNAL CBRT
68      INTEGER ig,l
69
70      LOGICAL firstcall,improved
71      SAVE firstcall
72
73      REAL zq(ngridmx,nlayermx,nqmx)
74      REAL zq0(ngridmx,nlayermx,nqmx)
75      REAL zqsat(ngridmx,nlayermx)
76      REAL zt(ngridmx,nlayermx)
77
78
79      REAL masse (ngridmx,nlayermx)
80      REAL epaisseur (ngridmx,nlayermx)
81      REAL dustcores(ngridmx,nlayermx)  !Dust number density (#/kg)
82      REAL rfinal        ! Ice crystal radius after condensation(m)
83      REAL seq           ! Equilibrium saturation ration (accounting for curvature effect)
84      REAL dzq           ! masse de glace echangee (kg/kg)
85      REAL lw       !Latent heat of sublimation (J.kg-1)
86      REAL To
87
88      REAL Ctot
89      REAL*8 ph2o,satu
90      REAL gr,Cste,up,dwn,newvap
91
92      DATA firstcall/.true./
93      DATA improved/.false./      ! Actionne une microphysique plus raffinee
94
95c     Reference temperature, T=273,15 K
96      data To/273.15/
97
98      SAVE improved
99
100c     POur diagnostique :
101c     ~~~~~~~~~~~~~~~~~
102      REAL taucond(ngridmx,nlayermx)   ! taux de condensation (kg/kg/s-1)
103      REAL mtot(ngridmx)               ! Total mass of water vapor (kg/m2)
104      REAL icetot(ngridmx)             ! Total mass of water ice (kg/m2)
105      REAL rave(ngridmx)               ! Mean crystal radius in a column (m)
106      REAL cloudmass(ngridmx,nlayermx) ! mass of ice in each layer (kg/m2)
107      REAL op825(ngridmx,nlayermx)    ! abs optical depth at 825 cm-1
108      REAL tau825(ngridmx)    ! abs optical depth at 825 cm-1
109      REAL waterhem(2,3)
110      real a,b, Qabs
111
112      REAL icetot2pm(ngridmx)
113      REAL mtot2pm(ngridmx)
114      REAL rave2pm(ngridmx)
115      REAL qsurf2pm(ngridmx)
116      INTEGER i,nit
117      REAL ecart,loctime
118      integer univtime,lon2pm(ngridmx)
119
120      integer ig_vl2
121      real latvl2,lonvl2
122
123      logical tesfield  ! output of TES like data
124      data tesfield /.false./
125
126c    ** un petit test de coherence
127c       --------------------------
128
129      IF (firstcall) THEN
130         IF(ngrid.NE.ngridmx) THEN
131            PRINT*,'STOP dans watercloud'
132            PRINT*,'probleme de dimensions :'
133            PRINT*,'ngrid  =',ngrid
134            PRINT*,'ngridmx  =',ngridmx
135            STOP
136         ENDIF
137         if(nq.gt.nqmx) then
138           write(*,*) 'stop in watercloud (nq.gt.nqmx)!'
139           write(*,*) 'nq=',nq,' nqmx=',nqmx
140           stop
141         endif
142         firstcall=.false.
143      ENDIF
144
145
146c-----------------------------------------------------------------------
147c    1. initialisation
148c    -----------------
149
150c    On "update" la valeur de q(nqmx) (water vapor) et temperature.
151c    On effectue qqes calculs preliminaires sur les couches :
152c    masse (kg.m-2), epaisseur(m).
153
154      do l=1,nlay
155          do ig=1,ngrid
156            zq(ig,l,nq)=pq(ig,l,nq)+pdq(ig,l,nq)*ptimestep
157            zq(ig,l,nq)=max(zq(ig,l,nq),1.E-30) ! FF 12/2004
158            zq0(ig,l,nq)=zq(ig,l,nq)
159            zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
160            masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
161            epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
162            if(iceparty) then
163              zq(ig,l,nq-1)=pq(ig,l,nq-1)+pdq(ig,l,nq-1)*ptimestep
164              zq(ig,l,nq-1)=max(zq(ig,l,nq-1),0.) ! FF 12/2004
165              zq0(ig,l,nq-1)=zq(ig,l,nq-1)
166              rdust(ig,l)= max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9)
167c         Calcul du rayon moyen des particules de glace.
168c         Hypothese : Dans une couche, la glace presente se
169c         repartit uniformement autour du nbre de poussieres
170c         dont le rayon moyen est prescrit par rdust.
171c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
172              if ((dustbin.eq.1).and..not.(doubleq)) then
173                zq(ig,l,1)=pq(ig,l,1)+pdq(ig,l,1)*ptimestep
174                dustcores(ig,l)=max(zq(ig,l,1)/
175     &           (rho_dust*4./3.*pi*radius(1)**3.),1.e-9)
176              else
177                dustcores(ig,l)=( epaisseur(ig,l)/masse(ig,l) ) *
178     &       2.e+6/0.1*max(tau(ig,1),0.001)*exp(-pzlay(ig,l)/10000.)
179
180c            TEMPORAIRE : réduction du nombre de nuclei FF 04/2008 :
181c               dustcores(ig,l) = dustcores(ig,l) / 27. ! reduction facteur 3
182                dustcores(ig,l) = dustcores(ig,l) / 8. ! reduction facteur 2
183              endif
184              rice(ig,l)=CBRT( ( zq(ig,l,nq-1)/rho_ice+
185     &         dustcores(ig,l)*(4./3.)*pi*rdust(ig,l)**3. )
186     &         / (dustcores(ig,l)*4./3.*pi) )
187              rice(ig,l)=max(rice(ig,l),rdust(ig,l))
188            endif
189          enddo
190      enddo
191
192      call zerophys(ngrid*nq,pdqscloud)
193      call zerophys(ngrid*nlay*nq,pdqcloud)
194      call zerophys(ngrid*nlay,pdtcloud)
195      call zerophys(ngrid,mtot)
196
197      call zerophys(ngrid,icetot)
198      call zerophys(ngrid,rave)
199      call zerophys(ngrid,cloudmass)
200
201c    ----------------------------------------------
202c
203c
204c       Rapport de melange a saturation dans la couche l : -------
205c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
206
207        call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
208
209c       taux de condensation (kg/kg/s-1) dans les differentes couches
210c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
211c       (Pour diagnostic seulement !)
212        if(.not.iceparty)then
213        do l=1, nlay
214          do ig=1,ngridmx
215            taucond(ig,l)= max((zq(ig,l,nq)-zqsat(ig,l))/ptimestep,0.)
216          end do
217        end do
218        endif
219
220
221        if(iceparty) then
222          do l=1,nlay
223            do ig=1,ngrid
224
225            If (improved) then
226c           Improved microphysics scheme
227c           ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
228
229             Ctot = zq(ig,l,nq) + zq(ig,l,nq-1)
230             ph2o = zq(ig,l,nq) * 44. / 18. * pplay(ig,l)
231             satu = zq(ig,l,nq) / zqsat(ig,l)
232
233             call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
234     &        ph2o,ph2o/satu,seq,rice(ig,l),gr)
235             Cste = ptimestep * 4. * pi * rice(ig,l)
236     *              * rho_ice * dustcores(ig,l)
237             up   = zq(ig,l,nq) + Cste * gr * seq
238             dwn  =    1.       + Cste * gr / zqsat(ig,l)
239             newvap = min(up/dwn,Ctot)
240
241             gr  = gr * ( newvap/zqsat(ig,l) - seq )
242             dzq = min( max( Cste * gr,-zq(ig,l,nq-1) )
243     *             , zq(ig,l,nq) )
244
245c            Nucleation (sat ratio must be larger than a critical value)
246c            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
247             if (satu.gt.1.) then
248               if (satu.le.1.4.and.zq(ig,l,nq-1).lt.1.e-8)
249     *            dzq = 0.
250             endif
251
252            Else
253c           Old version
254c           ~~~~~~~~~~~
255             if (zq(ig,l,nq).ge.zqsat(ig,l))then  !  Condensation
256               dzq=zq(ig,l,nq)-zqsat(ig,l)               
257             elseif(zq(ig,l,nq).lt.zqsat(ig,l))then  ! Sublimation
258               dzq=-min(zqsat(ig,l)-zq(ig,l,nq),zq(ig,l,nq-1))
259             endif
260
261            Endif
262
263c           Water Mass change
264c           ~~~~~~~~~~~~~~~~~
265            zq(ig,l,nq-1)=zq(ig,l,nq-1)+dzq
266            zq(ig,l,nq)=zq(ig,l,nq)-dzq
267
268            rice(ig,l)=max( CBRT ( (zq(ig,l,nq-1)/rho_ice
269     &       +dustcores(ig,l)*(4./3.)*pi*rdust(ig,l)**3.)
270     &       /(dustcores(ig,l)*4./3.*pi) ), rdust(ig,l))
271
272            if(activice.and.pclc(ig).gt.0)
273     &        rice(ig,l)=rice(ig,l)/CBRT(pclc(ig))
274             
275            enddo
276          enddo
277
278        else   ! if not iceparty
279
280c         Saturation couche nlay a 2 :
281c         ~~~~~~~~~~~~~~~~~~~~~~~~~~
282          do l=nlay,2, -1
283             do ig=1,ngrid
284             if (zq(ig,l,nq).gt.zqsat(ig,l))then
285               zq(ig,l-1,nq) =  zq(ig,l-1,nq)+(zq(ig,l,nq)-zqsat(ig,l))
286     &          *(pplev(ig,l)-pplev(ig,l+1))/(pplev(ig,l-1)-pplev(ig,l))
287               zq(ig,l,nq)=zqsat(ig,l)
288             endif
289             enddo
290          enddo
291
292c       Saturation couche l=1 si pas iceparty
293c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
294          do ig=1,ngridmx
295              if (zq(ig,1,nq).gt.zqsat(ig,1))then
296                 pdqscloud(ig,nq)= (zq(ig,1,nq)-zqsat(ig,1))
297     &         *(pplev(ig,1)-pplev(ig,2))/(g*ptimestep)
298                 zq(ig,1,nq)=zqsat(ig,1)
299              endif
300          end do
301
302        endif   ! (iceparty)
303
304c       Tendance finale
305c       ~~~~~~~~~~~~~~~
306        do l=1, nlay
307          do ig=1,ngridmx
308            pdqcloud(ig,l,nq) = (zq(ig,l,nq) - zq0(ig,l,nq))/ptimestep
309            if(iceparty) then
310              pdqcloud(ig,l,nq-1) =
311     &        (zq(ig,l,nq-1) - zq0(ig,l,nq-1))/ptimestep
312            endif
313            lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
314            pdtcloud(ig,l)=-pdqcloud(ig,l,nq)*lw/cpp
315          end do
316        end do
317
318c       A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
319c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
320c       Then that should not affect the ice particle radius
321        do ig=1,ngridmx
322          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
323            if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
324     &      rice(ig,2)=rice(ig,3)
325            rice(ig,1)=rice(ig,2)
326          end if
327        end do
328       
329c**************************************************
330c       Output 
331c**************************************************
332         do ig=1,ngridmx
333          do l=1 ,nlay
334c           masse de vapeur d'eau dans la couche l
335            mtot(ig)=mtot(ig)+pq(ig,l,nq)*masse(ig,l)
336            if(iceparty) then
337c             masse de glace d'eau dans la couche l
338              icetot(ig)=icetot(ig)+masse(ig,l)*pq(ig,l,nq-1)
339c             rayon moyen des cristaux dans la colonne ig
340              rave(ig)=rave(ig)+masse(ig,l)*pq(ig,l,nq-1)*rice(ig,l)
341              cloudmass(ig,l)=masse(ig,l)*pq(ig,l,nq-1)
342            endif
343          enddo
344          if(iceparty) then
345            rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
346            if (icetot(ig)*1000.lt.0.01) rave(ig)=0.
347          endif
348         enddo   ! (ngridmx)
349
350c         Computing abs optical depth at 825 cm-1 in each layer
351c         to simulate NEW TES retrieval
352          do ig=1,ngridmx
353            tau825(ig)=0.
354            do l=1 ,nlay
355               Qabs=min(max(0.4e6*rice(ig,l)-0.05 ,0.),1.2)
356               op825(ig,l)= 0.75*Qabs*pq(ig,l,nq-1)*masse(ig,l)
357     &                       /  (rho_ice*rice(ig,l))
358               tau825(ig)=tau825(ig)+ op825(ig,l)
359            end do
360          end do
361
362c ******************************************   
363c      Output in diagfi
364c ******************************************   
365
366      if (ngrid.gt.1) then
367
368          CALL WRITEDIAGFI(ngridmx,'mtot',
369     &  'total mass of water vapor','kg/m2',2,mtot)
370
371         if (callstats) then
372            call wstats(ngrid,"mtot",
373     .   "total mass of water vapor","kg/m2",2,mtot)
374         endif
375
376c        if(.not.iceparty)
377c     &  call WRITEDIAGFI(ngridmx,'taucond',
378c     &  'taux cond H2O ice','kg/kg/s',3,taucond)
379c
380         if(iceparty) then
381            CALL WRITEDIAGFI(ngridmx,'icetot',
382     &    'total mass of water ice','kg/m2',2,icetot)
383           if (callstats) then
384              call wstats(ngrid,"icetot",
385     .        "total mass of water ice","kg/m2",2,icetot)
386           endif
387c          CALL WRITEDIAGFI(ngridmx,'cloudmass',
388c     &    'mass of water ice/layer','kg/m2',3,cloudmass)
389c          CALL WRITEDIAGFI(ngridmx,'rice','ice radius',
390c     &    'meter',3,rice)
391           CALL WRITEDIAGFI(ngridmx,'rave','Mean ice radius',
392     &    'meter',2,rave)
393           CALL WRITEDIAGFI(ngridmx,'tauice','tau abs 825 cm-1',
394     &    '',2,tau825)
395         endif
396
397      else
398
399c          CALL WRITEG1D(ngridmx,1,mtot,'gas','kg/m2')
400c          CALL WRITEG1D(ngridmx,1,icetot,'ice','kg/m2')
401c          call WRITEG1D(ngridmx,nlay,rice,'rice','um')
402c          call WRITEG1D(ngridmx,nlay,dustcores,'pouss','#/layer')
403
404      endif
405 99   continue
406
407      RETURN
408      END
409 
Note: See TracBrowser for help on using the repository browser.