Context Navigation

watercloud_mod.F @ 2947

Last change on this file since 2947 was 2934, checked in by emillour, 21 months ago
Mars PCM: A first step in cleaning up outputs and adding field definitions for XIOS. EM
File size: 28.6 KB

Line
1	MODULE watercloud_mod
2
3	IMPLICIT NONE
4
5	REAL,SAVE,ALLOCATABLE :: zdqcloud(:,:,:) ! tendencies on pq due to condensation of H2O(kg/kg.s-1)
6	REAL,SAVE,ALLOCATABLE :: zdqscloud(:,:) ! tendencies on qsurf (calculated only by calchim but declared here)
7
8	!$OMP THREADPRIVATE(zdqcloud,zdqscloud)
9
10	CONTAINS
11
12	SUBROUTINE watercloud(ngrid,nlay,ptimestep,
13	& pplev,pplay,pdpsrf,pzlay,pt,pdt,
14	& pq,pdq,pdqcloud,pdtcloud,
15	& nq,tau,tauscaling,rdust,rice,nuice,
16	& rsedcloud,rhocloud,totcloudfrac)
17	USE ioipsl_getin_p_mod, ONLY : getin_p
18	USE updaterad, ONLY: updaterdust, updaterice_micro,
19	& updaterice_typ
20	USE improvedclouds_mod, ONLY: improvedclouds
21	USE watersat_mod, ONLY: watersat
22	use tracer_mod, only: nqmx, igcm_h2o_vap, igcm_h2o_ice,
23	& igcm_hdo_vap, igcm_hdo_ice,
24	& igcm_dust_mass, igcm_dust_number,
25	& igcm_ccn_mass, igcm_ccn_number,
26	& rho_dust, nuice_sed, nuice_ref,
27	& qperemin
28	use dimradmars_mod, only: naerkind
29	use write_output_mod, only: write_output
30	IMPLICIT NONE
31
32
33	c=======================================================================
34	c Water-ice cloud formation
35	c
36	c Includes two different schemes:
37	c - A simplified scheme (see simpleclouds.F)
38	c - An improved microphysical scheme (see improvedclouds.F)
39	c
40	c There is a time loop specific to cloud formation
41	c due to timescales smaller than the GCM integration timestep.
42	c
43	c Authors: Franck Montmessin, Francois Forget, Ehouarn Millour,
44	c J.-B. Madeleine, Thomas Navarro
45	c
46	c 2004 - 2012
47	c=======================================================================
48
49	c-----------------------------------------------------------------------
50	c declarations:
51	c -------------
52
53	include "callkeys.h"
54
55	c Inputs/outputs:
56	c ------
57
58	INTEGER, INTENT(IN) :: ngrid,nlay
59	INTEGER, INTENT(IN) :: nq ! nombre de traceurs
60	REAL, INTENT(IN) :: ptimestep ! pas de temps physique (s)
61	REAL, INTENT(IN) :: pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
62	REAL, INTENT(IN) :: pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
63	REAL, INTENT(IN) :: pdpsrf(ngrid) ! tendence surf pressure
64	REAL, INTENT(IN) :: pzlay(ngrid,nlay) ! altitude at the middle of the layers
65	REAL, INTENT(IN) :: pt(ngrid,nlay) ! temperature at the middle of the layers (K)
66	REAL, INTENT(IN) :: pdt(ngrid,nlay) ! tendence temperature des autres param.
67
68	REAL, INTENT(IN) :: pq(ngrid,nlay,nq) ! traceur (kg/kg)
69	rEAL, INTENT(IN) :: pdq(ngrid,nlay,nq) ! tendence avant condensation (kg/kg.s-1)
70
71	REAL, INTENT(IN) :: tau(ngrid,naerkind) ! Column dust optical depth at each point
72	REAL, INTENT(IN) :: tauscaling(ngrid) ! Convertion factor for dust amount
73	REAL, INTENT(INOUT) :: rdust(ngrid,nlay) ! Dust geometric mean radius (m)
74
75	REAL, INTENT(OUT) :: pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
76	REAL, INTENT(OUT) :: pdtcloud(ngrid,nlay) ! tendence temperature due
77	! a la chaleur latente
78	REAL, INTENT(INOUT) :: rice(ngrid,nlay) ! Ice mass mean radius (m)
79	! (r_c in montmessin_2004)
80	REAL, INTENT(OUT) :: nuice(ngrid,nlay) ! Estimated effective variance
81	! of the size distribution
82	REAL, INTENT(OUT) :: rsedcloud(ngrid,nlay) ! Cloud sedimentation radius
83	REAL, INTENT(OUT) :: rhocloud(ngrid,nlay) ! Cloud density (kg.m-3)
84
85	REAL, INTENT(INOUT):: totcloudfrac(ngrid) ! Cloud fraction (A. Pottier 2013)
86
87	c Locals:
88	c ------
89
90	! for ice radius computation
91	REAL Mo,No
92	REAl ccntyp
93
94	! for time loop
95	INTEGER microstep ! time subsampling step variable
96	INTEGER,SAVE :: imicro ! time subsampling for coupled water microphysics & sedimentation
97	REAL,SAVE :: microtimestep ! integration timestep for coupled water microphysics & sedimentation
98	REAL,SAVE :: microtimestep_prev=-999
99
100	!$OMP THREADPRIVATE(imicro,microtimestep)
101	!$OMP THREADPRIVATE(microtimestep_prev)
102	! tendency given by clouds (inside the micro loop)
103	REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud
104	REAL subpdtcloud(ngrid,nlay) ! cf. pdtcloud
105
106	! global tendency (clouds+physics)
107	REAL sum_subpdq(ngrid,nlay,nq) ! cf. pdqcloud
108	REAL sum_subpdt(ngrid,nlay) ! cf. pdtcloud
109
110	! no supersaturation when option supersat is false
111	REAL zt(ngrid,nlay) ! local value of temperature
112	REAL zqsat(ngrid,nlay) ! saturation
113
114	INTEGER iq,ig,l
115	LOGICAL,SAVE :: firstcall=.true.
116
117	!$OMP THREADPRIVATE(firstcall)
118	! HDO cycle
119	REAL :: alpha(ngrid,nlay) ! fractionation coefficient for HDO
120	REAL :: zq0(ngrid,nlay,nq) ! Initial mixing ratio: intermediate variable for HDO
121	! Representation of sub-grid water ice clouds A. Pottier 2013
122	REAL :: ztclf(ngrid, nlay)
123	REAL :: zqclf(ngrid, nlay,nq)
124	REAL :: zdelt
125	REAL :: norm
126	REAL :: ponder
127	REAL :: tcond(ngrid,nlay)
128	REAL :: zqvap(ngrid,nlay)
129	REAL :: zqice(ngrid,nlay)
130	REAL :: spant ! delta T for the temperature distribution
131	! REAL :: zqsat(ngrid,nlay) ! saturation
132	REAL :: pteff(ngrid, nlay)! effective temperature in the cloud,neb
133	REAL :: pqeff(ngrid, nlay, nq)! effective tracers quantities in the cloud
134	REAL :: cloudfrac(ngrid,nlay) ! cloud fraction
135	REAL :: mincloud ! min cloud frac
136	LOGICAL, save :: flagcloud=.true.
137
138	!$OMP THREADPRIVATE(flagcloud)
139
140	c ** un petit test de coherence
141	c --------------------------
142
143	IF (firstcall) THEN
144
145	if (nq.gt.nqmx) then
146	write(,) 'stop in watercloud (nq.gt.nqmx)!'
147	write(,) 'nq=',nq,' nqmx=',nqmx
148	call abort_physic("watercloud","nq.gt.nqmx",1)
149	endif
150
151	write(,) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
152	write(,) " igcm_h2o_ice=",igcm_h2o_ice
153
154	write(,) "time subsampling for microphysic ?"
155	#ifdef MESOSCALE
156	imicro = 2
157	#else
158	imicro = 30
159	#endif
160	call getin_p("imicro",imicro)
161	write(,)"watercloud: imicro = ",imicro
162
163	firstcall=.false.
164	ENDIF ! of IF (firstcall)
165
166	!! AS: moved out of firstcall to allow nesting+evoluting timestep
167	!! TBD: consider possible diff imicro with domains?
168	microtimestep = ptimestep/real(imicro)
169	if (microtimestep/=microtimestep_prev) then
170	! only tell the world if microtimestep has changed
171	write(,)"watercloud: Physical timestep is ",ptimestep
172	write(,)"watercloud: Microphysics timestep is ",microtimestep
173	microtimestep_prev=microtimestep
174	endif
175
176	c-----Initialization
177	sum_subpdq(1:ngrid,1:nlay,1:nq) = 0
178	sum_subpdt(1:ngrid,1:nlay) = 0
179
180	! default value if no ice
181	rhocloud(1:ngrid,1:nlay) = rho_dust
182
183	c-------------------------------------------------------------------
184	c 0. Representation of sub-grid water ice clouds
185	c------------------
186	c-----Initialization
187	pteff(1:ngrid,1:nlay) = 0
188	pqeff(1:ngrid,1:nlay,1:nq) = 0
189	DO l=1,nlay
190	DO ig=1,ngrid
191	pteff(ig,l)=pt(ig,l)
192	END DO
193	END DO
194	DO l=1,nlay
195	DO ig=1,ngrid
196	DO iq=1,nq
197	pqeff(ig,l,iq)=pq(ig,l,iq)
198	ENDDO
199	ENDDO
200	ENDDO
201	c-----Tendencies
202	DO l=1,nlay
203	DO ig=1,ngrid
204	ztclf(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
205	ENDDO
206	ENDDO
207	DO l=1,nlay
208	DO ig=1,ngrid
209	DO iq=1,nq
210	zqclf(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
211	ENDDO
212	ENDDO
213	ENDDO
214	c-----Effective temperature calculation
215	IF (CLFvarying) THEN
216	spant=3.0 ! delta T for the temprature distribution
217	mincloud=0.1 ! min cloudfrac when there is ice
218	IF (flagcloud) THEN
219	write(,) "Delta T", spant
220	write(,) "mincloud", mincloud
221	flagcloud=.false.
222	END IF
223	!CALL watersat(ngrid*nlay,ztclf,pplay,zqsat) !MV17: we dont need zqsat in the CLFvarying scheme
224	zqvap=zqclf(:,:,igcm_h2o_vap)
225	zqice=zqclf(:,:,igcm_h2o_ice)
226	CALL tcondwater(ngrid*nlay,pplay,zqvap+zqice,tcond)
227	DO l=1,nlay
228	DO ig=1,ngrid
229	zdelt=spant !MAX(spant*ztclf(ig,l),1.e-12), now totally in K. Fixed width
230	IF (tcond(ig,l) .ge. (ztclf(ig,l)+zdelt)) THEN
231	pteff(ig,l)=ztclf(ig,l)
232	cloudfrac(ig,l)=1.
233	ELSE IF (tcond(ig,l) .le. (ztclf(ig,l)-zdelt)) THEN
234	pteff(ig,l)=ztclf(ig,l)-zdelt
235	cloudfrac(ig,l)=mincloud
236	ELSE
237	cloudfrac(ig,l)=(tcond(ig,l)-ztclf(ig,l)+zdelt)/
238	& (2.0*zdelt)
239	pteff(ig,l)=(tcond(ig,l)+ztclf(ig,l)-zdelt)/2.
240	END IF
241	pteff(ig,l)=pteff(ig,l)-pdt(ig,l)*ptimestep
242	IF (cloudfrac(ig,l).le.mincloud) THEN !MV17: replaced .le.0 by .le.mincloud
243	cloudfrac(ig,l)=mincloud
244	ELSE IF (cloudfrac(ig,l).gt.1) THEN
245	cloudfrac(ig,l)=1.
246	END IF
247	ENDDO
248	ENDDO
249	c-----Calculation of the total cloud coverage of the column
250	DO ig=1,ngrid
251	totcloudfrac(ig) = 0.
252	norm=0.
253	DO l=1,nlay
254	ponder=zqice(ig,l)*(pplev(ig,l) - pplev(ig,l+1))
255	totcloudfrac(ig) = totcloudfrac(ig)
256	& + cloudfrac(ig,l)*ponder
257	norm=norm+ponder
258	ENDDO
259	totcloudfrac(ig)=MAX(totcloudfrac(ig)/norm,1.e-12) ! min value if NaNs
260	ENDDO
261	c-----Effective tracers quantities in the cloud fraction
262	IF (microphys) THEN
263	pqeff(:,:,igcm_ccn_mass)=pq(:,:,igcm_ccn_mass)/
264	& cloudfrac(:,:)
265	pqeff(:,:,igcm_ccn_number)=pq(:,:,igcm_ccn_number)/
266	& cloudfrac(:,:)
267	END IF ! end if (microphys)
268	pqeff(:,:,igcm_h2o_ice)=pq(:,:,igcm_h2o_ice)/
269	& cloudfrac(:,:)
270	!! CLFvarying outputs
271	! CALL write_output('pqeffice','pqeffice',
272	! & 'kg/kg',pqeff(:,:,igcm_h2o_ice))
273	! CALL write_output('pteff','pteff',
274	! & 'K',pteff(:,:))
275	! CALL write_output('tcond','tcond',
276	! & 'K',tcond(:,:))
277	! CALL write_output('cloudfrac','cloudfrac',
278	! & 'K',cloudfrac(:,:))
279	END IF ! end if (CLFvarying)
280	c------------------------------------------------------------------
281	c Time subsampling for microphysics
282	c------------------------------------------------------------------
283	rhocloud(1:ngrid,1:nlay) = rho_dust
284	DO microstep=1,imicro
285
286	c-------------------------------------------------------------------
287	c 1. Tendencies:
288	c------------------
289
290
291	c------ Temperature tendency subpdt
292	! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
293	! If imicro=1 subpdt is the same as pdt
294	DO l=1,nlay
295	DO ig=1,ngrid
296	sum_subpdt(ig,l) = sum_subpdt(ig,l)
297	& + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
298	ENDDO
299	ENDDO
300	c------ Tracers tendencies subpdq are additionned
301	c------ At each micro timestep we add pdq in order to have a stepped entry
302	IF (microphys) THEN
303	DO l=1,nlay
304	DO ig=1,ngrid
305	sum_subpdq(ig,l,igcm_dust_mass) =
306	& sum_subpdq(ig,l,igcm_dust_mass)
307	& + pdq(ig,l,igcm_dust_mass)
308	sum_subpdq(ig,l,igcm_dust_number) =
309	& sum_subpdq(ig,l,igcm_dust_number)
310	& + pdq(ig,l,igcm_dust_number)
311	sum_subpdq(ig,l,igcm_ccn_mass) =
312	& sum_subpdq(ig,l,igcm_ccn_mass)
313	& + pdq(ig,l,igcm_ccn_mass)
314	sum_subpdq(ig,l,igcm_ccn_number) =
315	& sum_subpdq(ig,l,igcm_ccn_number)
316	& + pdq(ig,l,igcm_ccn_number)
317	ENDDO
318	ENDDO
319	ENDIF
320	DO l=1,nlay
321	DO ig=1,ngrid
322	sum_subpdq(ig,l,igcm_h2o_ice) =
323	& sum_subpdq(ig,l,igcm_h2o_ice)
324	& + pdq(ig,l,igcm_h2o_ice)
325	sum_subpdq(ig,l,igcm_h2o_vap) =
326	& sum_subpdq(ig,l,igcm_h2o_vap)
327	& + pdq(ig,l,igcm_h2o_vap)
328	IF (hdo) THEN
329	sum_subpdq(ig,l,igcm_hdo_ice) =
330	& sum_subpdq(ig,l,igcm_hdo_ice)
331	& + pdq(ig,l,igcm_hdo_ice)
332	sum_subpdq(ig,l,igcm_hdo_vap) =
333	& sum_subpdq(ig,l,igcm_hdo_vap)
334	& + pdq(ig,l,igcm_hdo_vap)
335	ENDIF !hdo
336	ENDDO
337	ENDDO
338
339	c-------------------------------------------------------------------
340	c 2. Main call to the different cloud schemes:
341	c------------------------------------------------
342	IF (microphys) THEN
343	CALL improvedclouds(ngrid,nlay,microtimestep,
344	& pplay,pteff,sum_subpdt,
345	& pqeff,sum_subpdq,subpdqcloud,subpdtcloud,
346	& nq,tauscaling)
347
348	ELSE
349	CALL simpleclouds(ngrid,nlay,microtimestep,
350	& pplay,pzlay,pteff,sum_subpdt,
351	& pqeff,sum_subpdq,subpdqcloud,subpdtcloud,
352	& nq,tau,rice)
353	ENDIF
354
355	c-------------------------------------------------------------------
356	c 3. Updating tendencies after cloud scheme:
357	c-----------------------------------------------
358
359	IF (microphys) THEN
360	DO l=1,nlay
361	DO ig=1,ngrid
362	sum_subpdq(ig,l,igcm_dust_mass) =
363	& sum_subpdq(ig,l,igcm_dust_mass)
364	& + subpdqcloud(ig,l,igcm_dust_mass)
365	sum_subpdq(ig,l,igcm_dust_number) =
366	& sum_subpdq(ig,l,igcm_dust_number)
367	& + subpdqcloud(ig,l,igcm_dust_number)
368	sum_subpdq(ig,l,igcm_ccn_mass) =
369	& sum_subpdq(ig,l,igcm_ccn_mass)
370	& + subpdqcloud(ig,l,igcm_ccn_mass)
371	sum_subpdq(ig,l,igcm_ccn_number) =
372	& sum_subpdq(ig,l,igcm_ccn_number)
373	& + subpdqcloud(ig,l,igcm_ccn_number)
374	ENDDO
375	ENDDO
376	ENDIF
377	DO l=1,nlay
378	DO ig=1,ngrid
379	sum_subpdq(ig,l,igcm_h2o_ice) =
380	& sum_subpdq(ig,l,igcm_h2o_ice)
381	& + subpdqcloud(ig,l,igcm_h2o_ice)
382	sum_subpdq(ig,l,igcm_h2o_vap) =
383	& sum_subpdq(ig,l,igcm_h2o_vap)
384	& + subpdqcloud(ig,l,igcm_h2o_vap)
385
386	IF (hdo) THEN
387	sum_subpdq(ig,l,igcm_hdo_ice) =
388	& sum_subpdq(ig,l,igcm_hdo_ice)
389	& + subpdqcloud(ig,l,igcm_hdo_ice)
390	sum_subpdq(ig,l,igcm_hdo_vap) =
391	& sum_subpdq(ig,l,igcm_hdo_vap)
392	& + subpdqcloud(ig,l,igcm_hdo_vap)
393	ENDIF ! hdo
394
395	ENDDO
396	ENDDO
397
398	IF (activice) THEN
399	DO l=1,nlay
400	DO ig=1,ngrid
401	sum_subpdt(ig,l) =
402	& sum_subpdt(ig,l) + subpdtcloud(ig,l)
403	ENDDO
404	ENDDO
405	ENDIF
406
407	! !! Example of how to use writediagmicrofi useful to
408	! !! get outputs at each microphysical sub-timestep (better to be used in 1D)
409	! CALL WRITEDIAGMICROFI(ngrid,imicro,microstep,
410	! & microtimestep,'subpdtcloud',
411	! & 'subpdtcloud','K/s',1,subpdtcloud(:,:))
412
413	ENDDO ! of DO microstep=1,imicro
414
415	c-------------------------------------------------------------------
416	c 6. Compute final tendencies after time loop:
417	c------------------------------------------------
418
419	c------ Temperature tendency pdtcloud
420	DO l=1,nlay
421	DO ig=1,ngrid
422	pdtcloud(ig,l) =
423	& sum_subpdt(ig,l)/real(imicro)-pdt(ig,l)
424	ENDDO
425	ENDDO
426
427	c------ Tracers tendencies pdqcloud
428	DO l=1,nlay
429	DO ig=1,ngrid
430	pdqcloud(ig,l,igcm_h2o_ice) =
431	& sum_subpdq(ig,l,igcm_h2o_ice)/real(imicro)
432	& - pdq(ig,l,igcm_h2o_ice)
433	pdqcloud(ig,l,igcm_h2o_vap) =
434	& sum_subpdq(ig,l,igcm_h2o_vap)/real(imicro)
435	& - pdq(ig,l,igcm_h2o_vap)
436	IF (hdo) THEN
437	pdqcloud(ig,l,igcm_hdo_ice) =
438	& sum_subpdq(ig,l,igcm_hdo_ice)/real(imicro)
439	& - pdq(ig,l,igcm_hdo_ice)
440	pdqcloud(ig,l,igcm_hdo_vap) =
441	& sum_subpdq(ig,l,igcm_hdo_vap)/real(imicro)
442	& - pdq(ig,l,igcm_hdo_vap)
443	ENDIF !hdo
444	ENDDO
445	ENDDO
446
447	IF(microphys) THEN
448	DO l=1,nlay
449	DO ig=1,ngrid
450	pdqcloud(ig,l,igcm_ccn_mass) =
451	& sum_subpdq(ig,l,igcm_ccn_mass)/real(imicro)
452	& - pdq(ig,l,igcm_ccn_mass)
453	pdqcloud(ig,l,igcm_ccn_number) =
454	& sum_subpdq(ig,l,igcm_ccn_number)/real(imicro)
455	& - pdq(ig,l,igcm_ccn_number)
456	ENDDO
457	ENDDO
458	ENDIF
459
460	IF(scavenging) THEN
461	DO l=1,nlay
462	DO ig=1,ngrid
463	pdqcloud(ig,l,igcm_dust_mass) =
464	& sum_subpdq(ig,l,igcm_dust_mass)/real(imicro)
465	& - pdq(ig,l,igcm_dust_mass)
466	pdqcloud(ig,l,igcm_dust_number) =
467	& sum_subpdq(ig,l,igcm_dust_number)/real(imicro)
468	& - pdq(ig,l,igcm_dust_number)
469	ENDDO
470	ENDDO
471	ENDIF
472
473	c------- Due to stepped entry, other processes tendencies can add up to negative values
474	c------- Therefore, enforce positive values and conserve mass
475	IF(microphys) THEN
476	DO l=1,nlay
477	DO ig=1,ngrid
478	IF ((pq(ig,l,igcm_ccn_number) +
479	& ptimestep* (pdq(ig,l,igcm_ccn_number) +
480	& pdqcloud(ig,l,igcm_ccn_number)) .le. 1.)
481	& .or. (pq(ig,l,igcm_ccn_mass) +
482	& ptimestep* (pdq(ig,l,igcm_ccn_mass) +
483	& pdqcloud(ig,l,igcm_ccn_mass)) .le. 1.e-20)) THEN
484	pdqcloud(ig,l,igcm_ccn_number) =
485	& - pq(ig,l,igcm_ccn_number)/ptimestep
486	& - pdq(ig,l,igcm_ccn_number) + 1.
487	pdqcloud(ig,l,igcm_dust_number) =
488	& -pdqcloud(ig,l,igcm_ccn_number)
489	pdqcloud(ig,l,igcm_ccn_mass) =
490	& - pq(ig,l,igcm_ccn_mass)/ptimestep
491	& - pdq(ig,l,igcm_ccn_mass) + 1.e-20
492	pdqcloud(ig,l,igcm_dust_mass) =
493	& -pdqcloud(ig,l,igcm_ccn_mass)
494	ENDIF
495	ENDDO
496	ENDDO
497	ENDIF
498
499	IF(scavenging) THEN
500	DO l=1,nlay
501	DO ig=1,ngrid
502	IF ((pq(ig,l,igcm_dust_number) +
503	& ptimestep* (pdq(ig,l,igcm_dust_number) +
504	& pdqcloud(ig,l,igcm_dust_number)) .le. 1.)
505	& .or. (pq(ig,l,igcm_dust_mass) +
506	& ptimestep* (pdq(ig,l,igcm_dust_mass) +
507	& pdqcloud(ig,l,igcm_dust_mass)) .le. 1.e-20)) THEN
508	pdqcloud(ig,l,igcm_dust_number) =
509	& - pq(ig,l,igcm_dust_number)/ptimestep
510	& - pdq(ig,l,igcm_dust_number) + 1.
511	pdqcloud(ig,l,igcm_ccn_number) =
512	& -pdqcloud(ig,l,igcm_dust_number)
513	pdqcloud(ig,l,igcm_dust_mass) =
514	& - pq(ig,l,igcm_dust_mass)/ptimestep
515	& - pdq(ig,l,igcm_dust_mass) + 1.e-20
516	pdqcloud(ig,l,igcm_ccn_mass) =
517	& -pdqcloud(ig,l,igcm_dust_mass)
518	ENDIF
519	ENDDO
520	ENDDO
521	ENDIF
522
523	DO l=1,nlay
524	DO ig=1,ngrid
525
526	IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
527	& (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice))
528	& .le. qperemin) THEN
529	pdqcloud(ig,l,igcm_h2o_ice) =
530	& - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
531	pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
532	if (hdo) then
533	pdqcloud(ig,l,igcm_hdo_ice) =
534	& - pq(ig,l,igcm_hdo_ice)/ptimestep - pdq(ig,l,igcm_hdo_ice)
535	pdqcloud(ig,l,igcm_hdo_vap) = -pdqcloud(ig,l,igcm_hdo_ice)
536	endif
537	ENDIF
538	IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
539	& (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
540	& .le. qperemin) THEN
541	pdqcloud(ig,l,igcm_h2o_vap) =
542	& - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
543	pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
544	if (hdo) then
545	pdqcloud(ig,l,igcm_hdo_vap) =
546	& - pq(ig,l,igcm_hdo_vap)/ptimestep - pdq(ig,l,igcm_hdo_vap)
547	pdqcloud(ig,l,igcm_hdo_ice) = -pdqcloud(ig,l,igcm_hdo_vap)
548	endif
549	ENDIF
550
551	ENDDO
552	ENDDO
553
554	c------Update the ice and dust particle size "rice" for output or photochemistry
555	c------Only rsedcloud is used for the water cycle
556
557	IF(scavenging) THEN
558	DO l=1, nlay
559	DO ig=1,ngrid
560
561	call updaterdust(
562	& pq(ig,l,igcm_dust_mass) + ! dust mass
563	& (pdq(ig,l,igcm_dust_mass) + ! dust mass
564	& pdqcloud(ig,l,igcm_dust_mass))*ptimestep, ! dust mass
565	& pq(ig,l,igcm_dust_number) + ! dust number
566	& (pdq(ig,l,igcm_dust_number) + ! dust number
567	& pdqcloud(ig,l,igcm_dust_number))*ptimestep, ! dust number
568	& rdust(ig,l))
569
570	ENDDO
571	ENDDO
572	ENDIF
573
574	IF(microphys) THEN
575
576	! In case one does not want to allow supersatured water when using microphysics.
577	! Not done by default.
578	IF(.not.supersat) THEN
579	! !! initial mixing ratios for initial D/H ratio calculation
580	zq0(:,:,:) = pq(:,:,:) + pdq(:,:,:)*ptimestep
581	zt = pt + (pdt+pdtcloud)*ptimestep
582	call watersat(ngrid*nlay,zt,pplay,zqsat)
583	DO l=1, nlay
584	DO ig=1,ngrid
585	IF (pq(ig,l,igcm_h2o_vap)
586	& + (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
587	& * ptimestep .ge. zqsat(ig,l)) THEN
588	pdqcloud(ig,l,igcm_h2o_vap) =
589	& (zqsat(ig,l) - pq(ig,l,igcm_h2o_vap))/ptimestep
590	& - pdq(ig,l,igcm_h2o_vap)
591	pdqcloud(ig,l,igcm_h2o_ice) =
592	& -pdqcloud(ig,l,igcm_h2o_vap)
593	! !! HDO ice flux has to be modified in consequence
594	IF (hdo) THEN
595	! !! Logically only condensation can happen in this case
596	IF (pdqcloud(ig,l,igcm_h2o_ice) .gt. 0.0) THEN
597	IF ( zq0(ig,l,igcm_h2o_vap) .gt. qperemin ) THEN
598	! !! Lamb et al. 2017
599	alpha(ig,l) = exp(13525./zt(ig,l)**2.-5.59d-2)
600	pdqcloud(ig,l,igcm_hdo_ice) =
601	& pdqcloud(ig,l,igcm_h2o_ice)alpha(ig,l)
602	& ( zq0(ig,l,igcm_hdo_vap)
603	& /zq0(ig,l,igcm_h2o_vap) )
604	pdqcloud(ig,l,igcm_hdo_ice) =
605	& min(pdqcloud(ig,l,igcm_hdo_ice),
606	& zq0(ig,l,igcm_hdo_vap)/ptimestep)
607	ELSE
608	pdqcloud(ig,l,igcm_hdo_ice) = 0.
609	ENDIF
610	!! sublimation
611	ELSE
612	IF ( zq0(ig,l,igcm_h2o_ice) .gt. qperemin ) THEN
613	pdqcloud(ig,l,igcm_hdo_ice) =
614	& pdqcloud(ig,l,igcm_h2o_ice)*
615	& ( zq0(ig,l,igcm_hdo_ice)
616	& /zq0(ig,l,igcm_h2o_ice) )
617	pdqcloud(ig,l,igcm_hdo_ice) =
618	& max(pdqcloud(ig,l,igcm_hdo_ice),
619	& -zq0(ig,l,igcm_hdo_ice)/ptimestep)
620	ELSE
621	pdqcloud(ig,l,igcm_hdo_ice) = 0.
622	ENDIF
623	ENDIF !IF (pdqcloud(ig,l,igcm_h2o_ice).gt.0.)
624	pdqcloud(ig,l,igcm_hdo_vap) =
625	& -pdqcloud(ig,l,igcm_hdo_ice)
626	ENDIF !IF (hdo)
627	! no need to correct ccn_number, updaterad can handle this properly.
628	ENDIF
629	ENDDO
630	ENDDO
631	ENDIF
632
633	DO l=1, nlay
634	DO ig=1,ngrid
635
636	call updaterice_micro(
637	& pq(ig,l,igcm_h2o_ice) + ! ice mass
638	& (pdq(ig,l,igcm_h2o_ice) + ! ice mass
639	& pdqcloud(ig,l,igcm_h2o_ice))*ptimestep, ! ice mass
640	& pq(ig,l,igcm_ccn_mass) + ! ccn mass
641	& (pdq(ig,l,igcm_ccn_mass) + ! ccn mass
642	& pdqcloud(ig,l,igcm_ccn_mass))*ptimestep, ! ccn mass
643	& pq(ig,l,igcm_ccn_number) + ! ccn number
644	& (pdq(ig,l,igcm_ccn_number) + ! ccn number
645	& pdqcloud(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
646	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
647
648	ENDDO
649	ENDDO
650
651	ELSE ! no microphys
652
653	DO l=1,nlay
654	DO ig=1,ngrid
655
656	call updaterice_typ(
657	& pq(ig,l,igcm_h2o_ice) + ! ice mass
658	& (pdq(ig,l,igcm_h2o_ice) + ! ice mass
659	& pdqcloud(ig,l,igcm_h2o_ice))*ptimestep, ! ice mass
660	& tau(ig,1),pzlay(ig,l),rice(ig,l))
661
662	ENDDO
663	ENDDO
664
665	ENDIF ! of IF(microphys)
666
667
668
669	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
670	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
671	c Then that should not affect the ice particle radius
672	do ig=1,ngrid
673	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
674	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
675	& rice(ig,2)=rice(ig,3)
676	rice(ig,1)=rice(ig,2)
677	end if
678	end do
679
680
681	DO l=1,nlay
682	DO ig=1,ngrid
683	rsedcloud(ig,l)=max(rice(ig,l)*
684	& (1.+nuice_sed)(1.+nuice_sed)(1.+nuice_sed),
685	& rdust(ig,l))
686	! rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
687	ENDDO
688	ENDDO
689
690	! used for rad. transfer calculations
691	! nuice is constant because a lognormal distribution is prescribed
692	nuice(1:ngrid,1:nlay)=nuice_ref
693
694	c------Update tendencies for sub-grid water ice clouds
695	IF (CLFvarying) THEN
696	DO ig=1,ngrid
697	DO l=1,nlay
698	pdqcloud(ig,l,igcm_dust_mass)=pdqcloud(ig,l,igcm_dust_mass)
699	& *cloudfrac(ig,l)
700	pdqcloud(ig,l,igcm_ccn_mass)=pdqcloud(ig,l,igcm_ccn_mass)
701	& *cloudfrac(ig,l)
702	pdqcloud(ig,l,igcm_dust_number)=pdqcloud(ig,l,
703	& igcm_dust_number) *cloudfrac(ig,l)
704	pdqcloud(ig,l,igcm_ccn_number)=pdqcloud(ig,l,
705	& igcm_ccn_number) *cloudfrac(ig,l)
706	pdqcloud(ig,l,igcm_h2o_vap)=pdqcloud(ig,l,
707	& igcm_h2o_vap) *cloudfrac(ig,l)
708	pdqcloud(ig,l,igcm_h2o_ice)=pdqcloud(ig,l,
709	& igcm_h2o_ice) *cloudfrac(ig,l)
710	ENDDO
711	ENDDO
712	pdtcloud(:,:)=pdtcloud(:,:)*cloudfrac(:,:)
713	ENDIF
714	#ifndef MESOSCALE
715	c=======================================================================
716	call write_output("watercloud_pdqh2oice","pdqcloud_h2o_ice "//
717	& "after microphysics","kg/kg.s-1",pdqcloud(:,:,igcm_h2o_ice))
718	call write_output("watercloud_pdqh2ovap","pdqcloud_h2o_vap "//
719	& "after microphysics","kg/kg.s-1",pdqcloud(:,:,igcm_h2o_vap))
720	if (hdo) then
721	call write_output("watercloud_pdqhdoice","pdqcloud_hdo_ice "//
722	& "after microphysics","kg/kg.s-1",pdqcloud(:,:,igcm_hdo_ice))
723	call write_output("watercloud_pdqhdovap","pdqcloud_hdo_vap "//
724	& "after microphysics","kg/kg.s-1",pdqcloud(:,:,igcm_hdo_vap))
725	endif
726	! call write_output("pdqccn2","pdqcloudccn apres microphysique"
727	! & ,"kg/kg.s-1",pdqcloud(:,:,
728	! & igcm_ccn_mass))
729	! call write_output("pdqccnN2","pdqcloudccnN apres "//
730	! & "microphysique","nb/kg.s-1",pdqcloud(:,:,
731	! & igcm_ccn_number))
732	! call write_output("pdqdust2", "pdqclouddust apres "//
733	! & "microphysique","kg/kg.s-1",pdqcloud(:,:,
734	! & igcm_dust_mass))
735	! call write_output("pdqdustN2", "pdqclouddustN apres "//
736	! & "microphysique","nb/kg.s-1",pdqcloud(:,:,
737	! & igcm_dust_number))
738	c=======================================================================
739	#endif
740
741	END SUBROUTINE watercloud
742
743	subroutine ini_watercloud_mod(ngrid,nlayer,nq)
744	implicit none
745
746	integer,intent(in) :: ngrid ! number of atmospheric columns
747	integer,intent(in) :: nlayer ! number of atmospheric layers
748	integer,intent(in) :: nq ! number of tracers
749
750	allocate(zdqcloud(ngrid,nlayer,nq))
751	zdqcloud(:,:,:)=0
752	allocate(zdqscloud(ngrid,nq))
753	zdqscloud(:,:)=0
754
755	end subroutine ini_watercloud_mod
756
757
758	subroutine end_watercloud_mod
759	implicit none
760
761	if (allocated(zdqcloud)) deallocate(zdqcloud)
762	if (allocated(zdqscloud)) deallocate(zdqscloud)
763
764	end subroutine end_watercloud_mod
765
766	END MODULE watercloud_mod

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

Download in other formats:

Original Format