source: trunk/LMDZ.PLUTO.old/libf/phypluto/cocloud.F @ 3436

Last change on this file since 3436 was 3175, checked in by emillour, 11 months ago

Pluto PCM:
Add the old Pluto LMDZ for reference (required prior step to making
an LMDZ.PLUTO using the same framework as the other physics packages).
TB+EM
File size: 8.7 KB
RevLine 
[3175]1       SUBROUTINE cocloud(ngrid,nlay,naersize, ptimestep,
2     &                pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
3     &                pq,pdq,pdqcloud,pdqscloud,pdtcloud,
4     &                nq,rice_co)
5
6      use comgeomfi_h
7      IMPLICIT NONE
8
9c=======================================================================
10c     Treatment of saturation of CARBON MONOXIDE
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. CO condensee
16c    (si saturation) est remise dans la couche en dessous.
17c     CO condensee dans la couche du bas est deposee a la surface
18c       
19c
20c=======================================================================
21
22c-----------------------------------------------------------------------
23c   declarations:
24c   -------------
25
26#include "dimensions.h"
27#include "dimphys.h"
28#include "comcstfi.h"
29#include "callkeys.h"
30#include "tracer.h"
31
32c   Inputs:
33c   ------
34
35      INTEGER ngrid,nlay
36      REAL ptimestep             ! pas de temps physique (s)
37      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
38      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
39      REAL pdpsrf(ngrid)         ! tendance surf pressure
40      REAL pzlev(ngrid,nlay+1)   ! altitude at layer boundaries
41      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
42      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
43      REAL pdt(ngrid,nlay)       ! tendance temperature des autres param.
44
45      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
46      real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
47      integer naersize   ! nombre de traceurs radiativement actifs (=naerkind)
48      integer nq         ! nombre de traceurs
49
50c   Outputs:
51c   -------
52
53      real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation CO(kg/kg.s-1)
54      real pdqscloud(ngrid,nq)     ! flux en surface (kg.m-2.s-1)
55      REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
56                                   !   a la chaleur latente
57
58      REAL rice_co(ngrid,nlay)    ! Ice mass mean radius (m)
59                               ! (r_c in montmessin_2004)
60
61c   local:
62c   ------
63c      REAL Nmix   ! Cloud condensation nuclei
64c      parameter (Nmix=1.E2)  ! /kg
65!      parameter (Nmix=1)  ! /kg
66      real rnuclei  ! Nuclei geometric mean radius (m)
67      parameter (rnuclei=2.E-7)  ! m
68
69      REAL CBRT
70      EXTERNAL CBRT
71      INTEGER ig,l
72
73
74      REAL zq(ngridmx,nlayermx,nqmx)  ! local value of tracers
75      REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
76      REAL zqsat(ngridmx,nlayermx)    ! saturation
77      REAL zt(ngridmx,nlayermx)       ! local value of temperature
78
79      REAL vecnull(ngridmx*nlayermx)
80   
81      REAL masse (ngridmx,nlayermx)
82      REAL epaisseur (ngridmx,nlayermx)
83!      REAL rfinal        ! Ice crystal radius after condensation(m)
84      REAL*8 dzq           ! masse de glace echangee (kg/kg)
85      REAL lw       !Latent heat of sublimation (J.kg-1)
86
87
88      LOGICAL,SAVE :: firstcall=.true.
89
90! indexes of co gas, co ice and dust tracers:
91      INTEGER,SAVE :: i_co=0  ! co gas
92      INTEGER,SAVE :: i_ice=0  ! co ice
93
94c    ** un petit test de coherence
95c       --------------------------
96
97      IF (firstcall) THEN
98        IF(ngrid.NE.ngridmx) THEN
99            PRINT*,'STOP dans cocloud'
100            PRINT*,'probleme de dimensions :'
101            PRINT*,'ngrid  =',ngrid
102            PRINT*,'ngridmx  =',ngridmx
103            STOP
104        ENDIF
105         
106        if (nq.gt.nqmx) then
107           write(*,*) 'stop in cocloud (nq.gt.nqmx)!'
108           write(*,*) 'nq=',nq,' nqmx=',nqmx
109           stop
110        endif
111         
112c MELANIE : change these line
113        i_co=igcm_co_gas
114        i_ice=igcm_co_ice
115       
116        write(*,*) "cocloud: i_co=",i_co
117        write(*,*) "         i_ice=",i_ice
118
119        firstcall=.false.
120      ENDIF ! of IF (firstcall)
121
122
123c-----------------------------------------------------------------------
124c    1. initialisation
125c    -----------------
126
127c    On "update" la valeur de q(nqmx) (co vapor) et temperature.
128c    On effectue qqes calculs preliminaires sur les couches :
129c    masse (kg.m-2), epaisseur(m).
130
131      do l=1,nlay
132        do ig=1,ngrid
133          zq(ig,l,i_co)=pq(ig,l,i_co)+pdq(ig,l,i_co)*ptimestep
134          zq(ig,l,i_co)=max(zq(ig,l,i_co),1.E-30) ! FF 12/2004
135          zq0(ig,l,i_co)=zq(ig,l,i_co)
136          zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
137          masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
138          epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
139
140          zq(ig,l,i_ice)=pq(ig,l,i_ice)+pdq(ig,l,i_ice)*ptimestep
141          zq(ig,l,i_ice)=max(zq(ig,l,i_ice),0.) ! FF 12/2004
142          zq0(ig,l,i_ice)=zq(ig,l,i_ice)
143
144        enddo
145      enddo
146      pdqscloud(1:ngrid,1:nq)=0
147      pdqcloud(1:ngrid,1:nlay,1:nq)=0
148      pdtcloud(1:ngrid,1:nlay)=0
149
150c    ----------------------------------------------
151c
152c
153c       Rapport de melange a saturation dans la couche l : -------
154c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
155c         do l=1,nlay
156c          do ig=1,ngrid 
157c             call cosat(zt(ig,l),pplay(ig,l),zqsat(ig,l))
158c             write(101,*)'qsat',qsat(ig,l)
159c          enddo
160c         enddo
161
162         call cosat(ngridmx*nlayermx,zt,pplay,zqsat,vecnull,vecnull)
163c        TEMPORAIRE :
164c        test sans condensation atmospherique
165c        do l=1,nlay
166c         do ig=1,ngrid
167c           zqsat(ig,l) = zqsat(ig,l) *1000.
168c         end do
169c        end do
170c         call WRITEDIAGFI(ngrid,"qsat_co","qsat_co","unit",3,zqsat)
171c          do l=1,nlay
172c            do ig=1,ngrid
173c             zqsat(ig,l)=0.117*exp((16*568.7/8.314)*(1/90.7
174c     &         -1/zt(ig,l)))*100000
175c             zqsat(ig,l)=(zqsat(ig,l)/pplay(ig,l))*(16/28)
176c              write(106,*)'zqsat',zqsat(ig,l)
177c            enddo ! of do ig=1,ngrid
178c          enddo ! of do l=1,nlay
179
180c       taux de condensation (kg/kg/s-1) dans les differentes couches
181c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
182
183        do l=1,nlay
184          do ig=1,ngrid
185c           call cosat(zt(ig,l),pplay(ig,l),zqsat(ig,l))
186           if (zq(ig,l,i_co).ge.zqsat(ig,l))then  !  Condensation
187             dzq=zq(ig,l,i_co)-zqsat(ig,l)               
188           elseif(zq(ig,l,i_co).lt.zqsat(ig,l))then  ! Sublimation
189             dzq=-min(zqsat(ig,l)-zq(ig,l,i_co),zq(ig,l,i_ice))
190           endif
191
192c           CO Mass change
193c           ~~~~~~~~~~~~~~~~~
194            zq(ig,l,i_ice)=zq(ig,l,i_ice)+dzq
195            zq(ig,l,i_co)=zq(ig,l,i_co)-dzq
196            rice_co(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_co_ice
197     &       +Nmix_co*(4./3.)*pi*rnuclei**3.)
198     &       /(Nmix_co*4./3.*pi) ), rnuclei)  ! CBRT=cube root
199            enddo ! of do ig=1,ngrid
200          enddo ! of do l=1,nlay
201
202c         Saturation couche nlay a 2 :
203c         ~~~~~~~~~~~~~~~~~~~~~~~~~~
204c         do l=nlay,2, -1
205c          do ig=1,ngrid
206c           if (zq(ig,l,i_co).gt.zqsat(ig,l))then
207c             zq(ig,l-1,i_co)= zq(ig,l-1,i_co)+
208c     &                          (zq(ig,l,i_co)-zqsat(ig,l))
209c     &          *(pplev(ig,l)-pplev(ig,l+1))/(pplev(ig,l-1)-pplev(ig,l))
210c             zq(ig,l,i_co)=zqsat(ig,l)
211c           endif
212c          enddo
213c         enddo
214
215c       Saturation couche l=1 si pas iceparty
216c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
217c         do ig=1,ngridmx
218c           if (zq(ig,1,i_co).gt.zqsat(ig,1))then
219c             pdqscloud(ig,i_ice)=(zq(ig,1,i_co)-zqsat(ig,1))
220c     &           *(pplev(ig,1)-pplev(ig,2))/(g*ptimestep)
221c             zq(ig,1,i_co)=zqsat(ig,1)
222c           endif
223c         enddo
224
225c       Tendance finale
226c       ~~~~~~~~~~~~~~~
227        do l=1, nlay
228          do ig=1,ngridmx
229            pdqcloud(ig,l,i_co)=(zq(ig,l,i_co)
230     &                              -zq0(ig,l,i_co))/ptimestep
231            pdqcloud(ig,l,i_ice) =
232     &        (zq(ig,l,i_ice) - zq0(ig,l,i_ice))/ptimestep
233
234            lw=lw_co
235            pdtcloud(ig,l)=-pdqcloud(ig,l,i_co)*lw/cpp
236          end do
237        end do
238
239c       A correction if a lot of subliming co fills the 1st layer FF04/2005
240c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
241c       Then that should not affect the ice particle radius
242        do ig=1,ngridmx
243          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
244            if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
245     &      rice_co(ig,2)=rice_co(ig,3)
246            rice_co(ig,1)=rice_co(ig,2)
247          end if
248        end do
249       
250c**************************************************
251c       Output --- removed
252c**************************************************
253! NB: for diagnostics use zq(), the updated value of tracers
254c         Computing ext visible optical depth  in each layer
255
256      RETURN
257      END
258 
Note: See TracBrowser for help on using the repository browser.