source: trunk/LMDZ.PLUTO/libf/phypluto/ch4cloud.F90

Last change on this file was 3972, checked in by tbertrand, 5 weeks ago

PLUTO PCM:
Fixing bug in ch4clouds
TB

File size: 12.8 KB
Line 
1       SUBROUTINE ch4cloud(ngrid,nlay,ptimestep,  &
2                     pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt, &
3                     pq,pdq,pdqcloud,pdqscloud,pdtcloud,  &
4                     nq,rice_ch4)
5
6      use comgeomfi_h
7      use comcstfi_mod, only: pi, g, cpp, r
8      use tracer_h, only: igcm_ch4_gas, igcm_ch4_ice, &
9             rho_ch4_ice,lw_ch4,micro_indx
10      use callkeys_mod, only: Nmix_ch4, callmufi
11
12      IMPLICIT NONE
13
14!=======================================================================
15!     Treatment of saturation of METHANE
16!
17!
18!     Modif de zq si saturation dans l'atmosphere
19!     si zq(ig,l)> zqsat(ig,l) ->    zq(ig,l)=zqsat(ig,l)
20!     Le test est effectue de bas en haut. CH4 condensee
21!    (si saturation) est remise dans la couche en dessous.
22!     Le methane condensee dans la couche du bas est deposee a la surface
23!
24!     Melanie Vangvichith (adapted from Mars water clouds)
25!     Tanguy Bertrand
26!     Completely rewritten by Francois Forget (to properly estimate the
27!               latent heat release.
28!
29!=======================================================================
30
31!-----------------------------------------------------------------------
32!   declarations:
33!   -------------
34
35!   Inputs:
36!   ------
37
38      INTEGER ngrid,nlay
39      REAL ptimestep             ! pas de temps physique (s)
40      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
41      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
42      REAL pdpsrf(ngrid)         ! tendance surf pressure
43      REAL pzlev(ngrid,nlay+1)   ! altitude at layer boundaries
44      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
45      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
46      REAL pdt(ngrid,nlay)       ! tendance temperature des autres param.
47
48      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
49      real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
50      integer nq         ! nombre de traceurs
51
52!   Outputs:
53!   -------
54
55      real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation CH4(kg/kg.s-1)
56      real pdqscloud(ngrid,nq)     ! flux en surface (kg.m-2.s-1)
57      REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
58                                   !   a la chaleur latente
59      REAL rice_ch4(ngrid,nlay)    ! CH4 Ice mass mean radius (m)
60                               ! (r_c in montmessin_2004)
61
62!   local:
63!   ------
64!      REAL Nmix   ! Cloud condensation nuclei
65!      parameter (Nmix=1.E5)  ! /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(ngrid,nlay,nq)  ! local value of tracers
75      REAL zqsat(ngrid,nlay)    ! saturation
76      REAL zt(ngrid,nlay)       ! local value of temperature
77      !REAL rho(ngrid,nlay)       ! local value of atm density
78      REAL nmix(ngrid,nlay)     ! Profil of Nmix
79
80      REAL vecnull(ngrid*nlay)
81
82      LOGICAL,SAVE :: firstcall=.true.
83
84! indexes of methane gas, methane ice and dust tracers:
85      INTEGER,SAVE :: i_ch4=0  ! methane gas
86      INTEGER,SAVE :: i_ice=0  ! methane ice
87
88      REAL Tfin,qch4_fin
89      REAL temp_fin ! function used to compute Tfin at the end of the timestep
90
91!     TEMPORAIRE :
92
93      real  pdtcloudmax,pdtcloudmin
94      integer  igmin, igmax,lmin,lmax
95
96
97
98!    ** un petit test de coherence
99!       --------------------------
100
101      IF (firstcall) THEN
102        IF(ngrid.NE.ngrid) THEN
103            PRINT*,'STOP dans ch4cloud'
104            PRINT*,'probleme de dimensions :'
105            PRINT*,'ngrid  =',ngrid
106            PRINT*,'ngrid  =',ngrid
107            STOP
108        ENDIF
109
110        if (nq.gt.nq) then
111           write(*,*) 'stop in ch4cloud (nq.gt.nq)!'
112           write(*,*) 'nq=',nq,' nq=',nq
113           stop
114        endif
115
116        i_ch4=igcm_ch4_gas
117        i_ice=igcm_ch4_ice
118
119        if (i_ch4.eq.0) then
120          write(*,*) "stop in ch4cloud: i_ch4=0. "
121          write(*,*) "Maybe put ch4 in traceur.def?"
122          stop
123        endif
124
125        write(*,*) "methanecloud: i_ch4=",i_ch4
126        write(*,*) "            i_ice=",i_ice
127
128        firstcall=.false.
129      ENDIF ! of IF (firstcall)
130
131
132!-----------------------------------------------------------------------
133!    1. initialisation
134!    -----------------
135
136!    On "update" la valeur de q(nq) (methane vapor) et temperature.
137!    On effectue qqes calculs preliminaires sur les couches :
138
139      zt(:,:)=pt(:,:)+ pdt(:,:)*ptimestep
140      zq(:,:,i_ch4)=pq(:,:,i_ch4)+pdq(:,:,i_ch4)*ptimestep
141      zq(:,:,i_ice)=pq(:,:,i_ice)+pdq(:,:,i_ice)*ptimestep
142
143      !rho(:,:)=pplay(:,:)/(r*zt(:,:))
144      ! microphysics tracers are in particle/kg
145      if (callmufi) then
146        zq(:,:,micro_indx(1))=pq(:,:,micro_indx(1))+ &
147                    pdq(:,:,micro_indx(1))*ptimestep
148        zq(:,:,micro_indx(3))=pq(:,:,micro_indx(3))+ &
149                    pdq(:,:,micro_indx(3))*ptimestep
150      endif
151
152      do l=1,nlay
153        do ig=1,ngrid
154          zq(ig,l,i_ch4)=max(zq(ig,l,i_ch4),1.E-30)
155          zq(ig,l,i_ice)=max(zq(ig,l,i_ice),0.)
156        enddo
157      enddo
158      if (callmufi) then
159       do l=1,nlay
160        do ig=1,ngrid
161          zq(ig,l,micro_indx(1))=max(zq(ig,l,micro_indx(1)),1.E-30)
162          zq(ig,l,micro_indx(3))=max(zq(ig,l,micro_indx(3)),1.E-30)
163        enddo
164       enddo
165      endif
166
167      pdqscloud(1:ngrid,1:nq)=0.
168      pdqcloud(1:ngrid,1:nlay,1:nq)=0.
169      pdtcloud(1:ngrid,1:nlay)=0.
170
171      ! Nmix profile
172      if (callmufi) then
173        nmix(:,:)=zq(:,:,micro_indx(1))+zq(:,:,micro_indx(3))
174      else
175        nmix(:,:)=Nmix_ch4
176      endif
177!    ----------------------------------------------
178!
179!
180!       q à saturation dans la couche l à temperature prédite zt
181!       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
182         call methanesat(ngrid*nlay,zt,pplay,zqsat,vecnull)
183
184
185!       Loop to work where CH4 condensation/sublimation is occuring
186!       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
187        do l=1,nlay
188          do ig=1,ngrid
189             if ((zq(ig,l,i_ch4).gt.zqsat(ig,l)).or.  &
190            (zq(ig,l,i_ice).gt.1.E-10) ) then ! phase change
191
192               if(zq(ig,l,i_ice).lt.(zqsat(ig,l)-zq(ig,l,i_ch4))) then
193                 !case when everything is sublimed:
194                 qch4_fin = zq(ig,l,i_ch4) + zq(ig,l,i_ice)
195                 Tfin = zt(ig,l) - (qch4_fin-zq(ig,l,i_ch4))*lw_ch4/cpp
196               else
197                 !case when CH4 ice is present at the end
198                 qch4_fin = zqsat(ig,l)
199                 Tfin = zt(ig,l) - (qch4_fin-zq(ig,l,i_ch4))*lw_ch4/cpp
200                 if(abs(Tfin-zt(ig,l)).gt.1.) then
201!                  Improved calculation if the temperature change is significant (1K?)
202
203!                  Estimating the temperature Tfin and condensation at the end of the
204!                  timestep due to condensation-sublimation process, taking into acount
205!                  latent heat, given by:   Tfin - zt = (zq - qsat(Tfin))*L/cp
206!                  Tfin, solution, of the equation, is estimated by iteration
207!                  We assume that improved Tfin is between zt and the 1st estimation of Tfin:
208                   Tfin = temp_fin(zt(ig,l),Tfin,0.01,pplay(ig,l),  &
209                  zt(ig,l),zq(ig,l,i_ch4),qch4_fin)
210!                  Security check for small changes
211                   if(abs(zq(ig,l,i_ch4)-qch4_fin).gt.  &
212                     abs(zq(ig,l,i_ch4)- zqsat(ig,l))  ) then
213                      write(*,*) "warning: inaccuracy in CH4 clouds"
214                      write(*,*) 'zq(ig,l,i_ch4),qch4_fin, zqsat(ig,l)',  &
215                     zq(ig,l,i_ch4),qch4_fin, zqsat(ig,l)
216
217                      qch4_fin = zqsat(ig,l)
218
219                   end if
220                 end if
221               end if
222!              Final tendancies
223               pdqcloud(ig,l,i_ch4)=(qch4_fin-zq(ig,l,i_ch4))/ptimestep
224
225               !TB16
226               IF (zq(ig,l,i_ch4)+pdqcloud(ig,l,i_ch4)*ptimestep.lt.0.) &
227                                                      then
228                  pdqcloud(ig,l,i_ch4)=-zq(ig,l,i_ch4)/ptimestep
229               END IF
230
231               pdqcloud(ig,l,i_ice) = - pdqcloud(ig,l,i_ch4)
232               pdtcloud(ig,l)=(Tfin - zt(ig,l))/ptimestep
233
234!               Ice particle radius
235                zq(ig,l,i_ice)= &
236                 zq(ig,l,i_ice)+pdqcloud(ig,l,i_ice)*ptimestep
237                rice_ch4(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_ch4_ice &
238            +nmix(ig,l)*(4./3.)*pi*rnuclei**3.)/(nmix(ig,l)*4./3.*pi)), &
239                                             rnuclei)
240             end if
241            enddo ! of do ig=1,ngrid
242          enddo ! of do l=1,nlay
243
244!       A correction if a lot of subliming N2 fills the 1st layer FF04/2005
245!       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
246!       Then that should not affect the ice particle radius
247        do ig=1,ngrid
248          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
249            if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3))) &
250           rice_ch4(ig,2)=rice_ch4(ig,3)
251            rice_ch4(ig,1)=rice_ch4(ig,2)
252          end if
253        end do
254
255!************ ANALYSE pour Icarus ****************
256!
257!        pdtcloudmax =0.
258!        pdtcloudmin =0.
259!        igmin=999
260!        igmax=999
261!        lmin =999
262!        lmax =999
263!        do l=1, nlay
264!          do ig=1,ngrid
265!            if(pdtcloud(ig,l).gt.pdtcloudmax) then
266!               igmax=ig
267!               lmax = l
268!               pdtcloudmax=pdtcloud(ig,l)
269!            else if(pdtcloud(ig,l).lt.pdtcloudmin) then
270!               igmin=ig
271!               lmin = l
272!               pdtcloudmin=pdtcloud(ig,l)
273!            end if
274!          end do
275!        end do
276!
277!        write(*,*)
278!        write(*,*) "MAX , ig, l, Tbefore, Tafter, Tfin_ini"
279!        write(*,*) igmax,lmax, zt(igmax,lmax),
280!    &      zt(igmax,lmax)+ pdtcloud(igmax,lmax)*ptimestep,
281!    &zt(igmax,lmax)-(zqsat(igmax,lmax)-zq(igmax,lmax,i_ch4))*lw_ch4/cpp
282!        write(*,*) "MAX ,  Qch4gas_before, Qch4gas_after, Qsat ini"
283!       write(*,*) pq(igmax,lmax,i_ch4)+pdq(igmax,lmax,i_ch4)*ptimestep,
284!    &   pq(igmax,lmax,i_ch4)+
285!    &   (pdq(igmax,lmax,i_ch4)+ pdqcloud(igmax,lmax,i_ch4))*ptimestep
286!    &    ,zqsat(igmax,lmax)
287!        write(*,*) "MAX ,  Qch4ice_before, Qch4ice_after"
288!       write(*,*) pq(igmax,lmax,i_ice)+pdq(igmax,lmax,i_ice)*ptimestep,
289!    &   pq(igmax,lmax,i_ice)+
290!    &   (pdq(igmax,lmax,i_ice)+ pdqcloud(igmax,lmax,i_ice))*ptimestep
291!        write(*,*) 'Rice=',  rice_ch4(igmax,lmax)
292!
293!        write(*,*)
294!        write(*,*) "MIN , ig, l, Tbefore, Tafter, Tfin_ini"
295!        write(*,*) igmin,lmin, zt(igmin,lmin),
296!    &             zt(igmin,lmin)+ pdtcloud(igmin,lmin)*ptimestep,
297!    &zt(igmin,lmin)-(zqsat(igmin,lmin)-zq(igmin,lmin,i_ch4))*lw_ch4/cpp
298!        write(*,*) "MIN ,  Qch4gas_before, Qch4gas_after, Qsat ini"
299!       write(*,*) pq(igmin,lmin,i_ch4)+pdq(igmin,lmin,i_ch4)*ptimestep,
300!    &   pq(igmin,lmin,i_ch4)+
301!    &   (pdq(igmin,lmin,i_ch4)+ pdqcloud(igmin,lmin,i_ch4))*ptimestep
302!    &    ,zqsat(igmin,lmin)
303!        write(*,*) "MIN ,  Qch4ice_before, Qch4ice_after"
304!       write(*,*) pq(igmin,lmin,i_ice)+pdq(igmin,lmin,i_ice)*ptimestep,
305!    &   pq(igmin,lmin,i_ice)+
306!    &   (pdq(igmin,lmin,i_ice)+ pdqcloud(igmin,lmin,i_ice))*ptimestep
307!        write(*,*) 'Rice=',  rice_ch4(igmin,lmin)
308!        write(*,*) "pdtcloud(igmin,lmin) " , pdtcloud(igmin,lmin)
309!        write(*,*) "pdqcloud(i_ch4)", pdqcloud(igmin,lmin,i_ch4)
310!        write(*,*) "pdqcloud(i_ice)", pdqcloud(igmin,lmin,i_ice)
311!
312
313!************ FIN ANALYSE pour Icarus ****************
314
315!**************************************************
316!       Output --- removed
317!**************************************************
318! NB: for diagnostics use zq(), the updated value of tracers
319!         Computing ext visible optical depth  in each layer
320
321      RETURN
322      END
323
324!================================================================================
325
326      FUNCTION temp_fin(X1,X2,XACC,pres,temp,zq,qsat)
327       use comcstfi_mod, only: pi, g, cpp
328       use tracer_h, only: igcm_ch4_gas, igcm_ch4_ice, rho_ch4_ice, lw_ch4
329       implicit none
330
331       real pres, temp, zq,qsat
332       real X1,X2,XACC, F
333       real FMID,D,temp_fin,DX, XMID
334       integer JMAX,J
335       real func
336
337      PARAMETER (JMAX=40)
338
339      call methanesat(1,X2,pres,qsat,0.)
340      FMID = X2 -temp -(zq-qsat)*lw_ch4/cpp
341
342      call methanesat(1,X1,pres,qsat,0.)
343      F = X1 -temp -(zq-qsat)*lw_ch4/cpp
344
345      IF(F*FMID.GE.0.) write(*,*)  'Fix Tfin firt guesses in ch4cloud'
346      IF(F.LT.0.)THEN
347        temp_fin=X1
348        DX=X2-X1
349      ELSE
350        temp_fin=X2
351        DX=X1-X2
352      ENDIF
353      DO 11 J=1,JMAX
354        DX=DX*.5
355        XMID=temp_fin+DX
356
357        call methanesat(1,XMID,pres,qsat,0.)
358        FMID = XMID -temp -(zq-qsat)*lw_ch4/cpp
359
360        IF(FMID.LE.0.)temp_fin=XMID
361        IF(ABS(DX).LT.XACC .OR. FMID.EQ.0.) RETURN
36211    CONTINUE
363!     PAUSE 'too many bisections'
364      END
365
366
Note: See TracBrowser for help on using the repository browser.