Context Navigation

watercloud.F @ 635

Last change on this file since 635 was 633, checked in by tnavarro, 13 years ago
last scheme in commit r626 led to a wrong physical behaviour. This version uses a new subtimestep for microphysics that should be faster than the last one.
File size: 14.2 KB

Line
1	SUBROUTINE watercloud(ngrid,nlay,ptimestep,
2	& pplev,pplay,pdpsrf,pzlay,pt,pdt,
3	& pq,pdq,pdqcloud,pdtcloud,
4	& nq,tau,tauscaling,rdust,rice,nuice,
5	& rsedcloud,rhocloud)
6	! to use 'getin'
7	USE ioipsl_getincom
8	IMPLICIT NONE
9
10
11	c=======================================================================
12	c Water-ice cloud formation
13	c
14	c Includes two different schemes:
15	c - A simplified scheme (see simpleclouds.F)
16	c - An improved microphysical scheme (see improvedclouds.F)
17	c
18	c There is a time loop specific to cloud formation
19	c due to timescales smaller than the GCM integration timestep.
20	c
21	c Authors: Franck Montmessin, Francois Forget, Ehouarn Millour,
22	c J.-B. Madeleine, Thomas Navarro
23	c
24	c 2004 - 2012
25	c=======================================================================
26
27	c-----------------------------------------------------------------------
28	c declarations:
29	c -------------
30
31	#include "dimensions.h"
32	#include "dimphys.h"
33	#include "comcstfi.h"
34	#include "callkeys.h"
35	#include "tracer.h"
36	#include "comgeomfi.h"
37	#include "dimradmars.h"
38
39	c Inputs:
40	c ------
41
42	INTEGER ngrid,nlay
43	INTEGER nq ! nombre de traceurs
44	REAL ptimestep ! pas de temps physique (s)
45	REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
46	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
47	REAL pdpsrf(ngrid) ! tendence surf pressure
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
59	c Outputs:
60	c -------
61
62	real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
63	REAL pdtcloud(ngrid,nlay) ! tendence temperature due
64	! a la chaleur latente
65
66	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
67	! (r_c in montmessin_2004)
68	REAL nuice(ngrid,nlay) ! Estimated effective variance
69	! of the size distribution
70	real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius
71	real rhocloud(ngridmx,nlay) ! Cloud density (kg.m-3)
72
73	c local:
74	c ------
75
76	! for ice radius computation
77	REAL Mo,No
78	REAl ccntyp
79
80	! for time loop
81	INTEGER microstep ! time subsampling step variable
82	INTEGER imicro ! time subsampling for coupled water microphysics & sedimentation
83	SAVE imicro
84	REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
85	SAVE microtimestep
86
87	! tendency given by clouds (inside the micro loop)
88	REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud
89	REAL subpdtcloud(ngrid,nlay) ! cf. pdtcloud
90
91	! global tendency (clouds+physics)
92	REAL subpdq(ngrid,nlay,nq) ! cf. pdqcloud
93	REAL subpdt(ngrid,nlay) ! cf. pdtcloud
94
95	REAL CBRT
96	EXTERNAL CBRT
97
98
99	INTEGER iq,ig,l
100	LOGICAL,SAVE :: firstcall=.true.
101
102	c ** un petit test de coherence
103	c --------------------------
104
105	IF (firstcall) THEN
106	IF(ngrid.NE.ngridmx) THEN
107	PRINT*,'STOP dans watercloud'
108	PRINT*,'probleme de dimensions :'
109	PRINT*,'ngrid =',ngrid
110	PRINT*,'ngridmx =',ngridmx
111	STOP
112	ENDIF
113
114	if (nq.gt.nqmx) then
115	write(,) 'stop in watercloud (nq.gt.nqmx)!'
116	write(,) 'nq=',nq,' nqmx=',nqmx
117	stop
118	endif
119
120	write(,) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
121	write(,) " igcm_h2o_ice=",igcm_h2o_ice
122
123	write(,) "time subsampling for microphysic ?"
124	#ifdef MESOSCALE
125	imicro = 2
126	#else
127	imicro = 15
128	#endif
129	call getin("imicro",imicro)
130	write(,)"imicro = ",imicro
131
132	microtimestep = ptimestep/real(imicro)
133	write(,)"Physical timestep is",ptimestep
134	write(,)"Microphysics timestep is",microtimestep
135
136	firstcall=.false.
137	ENDIF ! of IF (firstcall)
138
139	c-----Initialization
140	subpdq(1:ngrid,1:nlay,1:nq) = 0
141	subpdt(1:ngrid,1:nlay) = 0
142
143	! default value if no ice
144	rhocloud(1:ngrid,1:nlay) = rho_dust
145
146
147
148	c------------------------------------------------------------------
149	c Time subsampling for microphysics
150	c------------------------------------------------------------------
151	DO microstep=1,imicro
152
153	c-------------------------------------------------------------------
154	c 1. Tendencies:
155	c------------------
156
157
158	c------ Temperature tendency subpdt
159	! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
160	! If imicro=1 subpdt is the same as pdt
161	DO l=1,nlay
162	DO ig=1,ngrid
163	subpdt(ig,l) = subpdt(ig,l)
164	& + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
165	ENDDO
166	ENDDO
167	c------ Tracers tendencies subpdq
168	c------ At each micro timestep we add pdq in order to have a stepped entry
169	IF (microphys) THEN
170	DO l=1,nlay
171	DO ig=1,ngrid
172	subpdq(ig,l,igcm_dust_mass) =
173	& subpdq(ig,l,igcm_dust_mass)
174	& + pdq(ig,l,igcm_dust_mass)
175	subpdq(ig,l,igcm_dust_number) =
176	& subpdq(ig,l,igcm_dust_number)
177	& + pdq(ig,l,igcm_dust_number)
178	subpdq(ig,l,igcm_ccn_mass) =
179	& subpdq(ig,l,igcm_ccn_mass)
180	& + pdq(ig,l,igcm_ccn_mass)
181	subpdq(ig,l,igcm_ccn_number) =
182	& subpdq(ig,l,igcm_ccn_number)
183	& + pdq(ig,l,igcm_ccn_number)
184	ENDDO
185	ENDDO
186	ENDIF
187	DO l=1,nlay
188	DO ig=1,ngrid
189	subpdq(ig,l,igcm_h2o_ice) =
190	& subpdq(ig,l,igcm_h2o_ice)
191	& + pdq(ig,l,igcm_h2o_ice)
192	subpdq(ig,l,igcm_h2o_vap) =
193	& subpdq(ig,l,igcm_h2o_vap)
194	& + pdq(ig,l,igcm_h2o_vap)
195	ENDDO
196	ENDDO
197
198
199	c-------------------------------------------------------------------
200	c 2. Main call to the different cloud schemes:
201	c------------------------------------------------
202	IF (microphys) THEN
203	CALL improvedclouds(ngrid,nlay,microtimestep,
204	& pplay,pt,subpdt,
205	& pq,subpdq,subpdqcloud,subpdtcloud,
206	& nq,tauscaling)
207
208	ELSE
209	CALL simpleclouds(ngrid,nlay,microtimestep,
210	& pplay,pzlay,pt,subpdt,
211	& pq,subpdq,subpdqcloud,subpdtcloud,
212	& nq,tau)
213	ENDIF
214
215
216	c-------------------------------------------------------------------
217	c 3. Updating tendencies after cloud scheme:
218	c-----------------------------------------------
219
220	IF (microphys) THEN
221	DO l=1,nlay
222	DO ig=1,ngrid
223	subpdq(ig,l,igcm_dust_mass) =
224	& subpdq(ig,l,igcm_dust_mass)
225	& + subpdqcloud(ig,l,igcm_dust_mass)
226	subpdq(ig,l,igcm_dust_number) =
227	& subpdq(ig,l,igcm_dust_number)
228	& + subpdqcloud(ig,l,igcm_dust_number)
229	subpdq(ig,l,igcm_ccn_mass) =
230	& subpdq(ig,l,igcm_ccn_mass)
231	& + subpdqcloud(ig,l,igcm_ccn_mass)
232	subpdq(ig,l,igcm_ccn_number) =
233	& subpdq(ig,l,igcm_ccn_number)
234	& + subpdqcloud(ig,l,igcm_ccn_number)
235	ENDDO
236	ENDDO
237	ENDIF
238	DO l=1,nlay
239	DO ig=1,ngrid
240	subpdq(ig,l,igcm_h2o_ice) =
241	& subpdq(ig,l,igcm_h2o_ice)
242	& + subpdqcloud(ig,l,igcm_h2o_ice)
243	subpdq(ig,l,igcm_h2o_vap) =
244	& subpdq(ig,l,igcm_h2o_vap)
245	& + subpdqcloud(ig,l,igcm_h2o_vap)
246	ENDDO
247	ENDDO
248
249
250	ENDDO ! of DO microstep=1,imicro
251
252	c-------------------------------------------------------------------
253	c 6. Compute final tendencies after time loop:
254	c------------------------------------------------
255
256	c------ Temperature tendency pdtcloud
257	DO l=1,nlay
258	DO ig=1,ngrid
259	pdtcloud(ig,l) =
260	& subpdt(ig,l)/real(imicro)-pdt(ig,l)
261	ENDDO
262	ENDDO
263	c------ Tracers tendencies pdqcloud
264	DO iq=1,nq
265	DO l=1,nlay
266	DO ig=1,ngrid
267	pdqcloud(ig,l,iq) = subpdq(ig,l,iq)/real(imicro)
268	& - pdq(ig,l,iq)
269	ENDDO
270	ENDDO
271	ENDDO
272
273	c------- Due to stepped entry, other processes tendencies can add up to negative values
274	c------- Therefore, enforce positive values and conserve mass
275
276	IF(microphys) THEN
277	DO l=1,nlay
278	DO ig=1,ngrid
279	IF (pq(ig,l,igcm_ccn_number) +
280	& ptimestep* (pdq(ig,l,igcm_ccn_number) +
281	& pdqcloud(ig,l,igcm_ccn_number)) .le. 0.) THEN
282	pdqcloud(ig,l,igcm_ccn_number) =
283	& - pq(ig,l,igcm_ccn_number)/ptimestep
284	& - pdq(ig,l,igcm_ccn_number)
285	pdqcloud(ig,l,igcm_dust_number) =
286	& -pdqcloud(ig,l,igcm_ccn_number)
287	pdqcloud(ig,l,igcm_ccn_mass) =
288	& - pq(ig,l,igcm_ccn_mass)/ptimestep
289	& - pdq(ig,l,igcm_ccn_mass)
290	pdqcloud(ig,l,igcm_dust_mass) =
291	& -pdqcloud(ig,l,igcm_ccn_mass)
292	ENDIF
293	ENDDO
294	ENDDO
295	ENDIF
296
297	DO l=1,nlay
298	DO ig=1,ngrid
299	IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
300	& (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice))
301	& .le. 1.e-8) THEN
302	pdqcloud(ig,l,igcm_h2o_ice) =
303	& - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
304	pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
305	ENDIF
306	IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
307	& (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
308	& .le. 1.e-8) THEN
309	pdqcloud(ig,l,igcm_h2o_vap) =
310	& - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
311	pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
312	ENDIF
313	ENDDO
314	ENDDO
315
316
317	c------Update the ice and dust particle size "rice" for output or photochemistry
318	c------Only rsedcloud is used for the water cycle
319
320	IF(scavenging) THEN
321	DO l=1, nlay
322	DO ig=1,ngrid
323
324	rdust(ig,l)=
325	& CBRT(r3n_q*(pq(ig,l,igcm_dust_mass)+
326	& (pdq(ig,l,igcm_dust_mass)
327	& +pdqcloud(ig,l,igcm_dust_mass))*ptimestep)/
328	& (pq(ig,l,igcm_dust_number)+
329	& (pdq(ig,l,igcm_dust_number)+
330	& pdqcloud(ig,l,igcm_dust_number)*ptimestep)))
331	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
332
333	Mo = pq(ig,l,igcm_h2o_ice)
334	& + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep
335	& + (pq(ig,l,igcm_ccn_mass)
336	& + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
337	& *tauscaling(ig) + 1.e-30
338	No = (pq(ig,l,igcm_ccn_number)
339	& + pdqcloud(ig,l,igcm_ccn_number)*ptimestep)
340	& *tauscaling(ig) + 1.e-30
341	rhocloud(ig,l) = (pq(ig,l,igcm_h2o_ice) +
342	& pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)
343	& / Mo * rho_ice
344	& + (pq(ig,l,igcm_ccn_mass)
345	& + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
346	& tauscaling(ig)/ Mo rho_dust
347	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
348	if ((Mo.lt.1.e-15) .or. (No.le.1.)) then
349	rice(ig,l) = 1.e-8
350	else
351	!! AS: only perform computations if conditions not met
352	rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.)
353	rice(ig,l)=min(max(rice(ig,l),1.e-10),500.e-6)
354	endif
355	rsedcloud(ig,l)=max( rice(ig,l)(1.+nuice_sed)
356	& (1.+nuice_sed)*(1.+nuice_sed),
357	& rdust(ig,l) )
358	rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
359	ENDDO
360	ENDDO
361	ELSE
362	DO l=1,nlay
363	DO ig=1,ngrid
364
365	rdust(ig,l)=
366	& CBRT(r3n_q*(pq(ig,l,igcm_dust_mass)+
367	& (pdq(ig,l,igcm_dust_mass)
368	& +pdqcloud(ig,l,igcm_dust_mass))*ptimestep)/
369	& (pq(ig,l,igcm_dust_number)+
370	& (pdq(ig,l,igcm_dust_number)+
371	& pdqcloud(ig,l,igcm_dust_number)*ptimestep)))
372	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
373
374	ccntyp =
375	& 1.3e+8max(tau(ig,1),0.001)/0.1exp(-pzlay(ig,l)/10000.)
376	ccntyp = ccntyp /ccn_factor
377	rice(ig,l)=max( CBRT ( ((pq(ig,l,igcm_h2o_ice)
378	& + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)/rho_ice
379	& +ccntyp(4./3.)pirdust(ig,l)rdust(ig,l)*rdust(ig,l))
380	& /(ccntyp4./3.pi) ), rdust(ig,l))
381	rsedcloud(ig,l)=max( rice(ig,l)(1.+nuice_sed)
382	& (1.+nuice_sed)*(1.+nuice_sed),
383	& rdust(ig,l) )
384	rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
385	ENDDO
386	ENDDO
387	ENDIF ! of IF(scavenging)
388
389
390	! used for rad. transfer calculations
391	! nuice is constant because a lognormal distribution is prescribed
392	nuice(1:ngrid,1:nlay)=nuice_ref
393
394
395	!--------------------------------------------------------------
396	!--------------------------------------------------------------
397
398
399	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
400	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
401	c Then that should not affect the ice particle radius
402	do ig=1,ngridmx
403	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
404	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
405	& rice(ig,2)=rice(ig,3)
406	rice(ig,1)=rice(ig,2)
407	end if
408	end do
409
410	c=======================================================================
411
412	END
413

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format