Context Navigation

watercloud_mod.F @ 2932

Last change on this file since 2932 was 2932, checked in by romain.vande, 21 months ago

Mars PCM:
Add a new routine to write output called write_output.
It needs to be imported (for example : use write_output_mod, only: write_output)
Then, it requires only 4 arguments : the name of the variable, its title, its units and the variable itself.
It detects the dimension of the variable and decide to output either in diagfi or diagsoil.
It is also compatible with XIOS (xml file needs to be adapted)
Writediagfi and writediagsoil routines are still available in case.
RV

File size: 28.5 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("pdqice2","pdqcloudice apres microphysique"
717	& ,"kg/kg.s-1",pdqcloud(:,:,igcm_h2o_ice))
718	call write_output("pdqvap2","pdqcloudvap apres microphysique"
719	& ,"kg/kg.s-1",pdqcloud(:,:,
720	& igcm_h2o_vap))
721	if (hdo) then
722	call write_output("pdqiceD","pdqiceD apres microphysique"
723	& ,"kg/kg.s-1",pdqcloud(:,:,igcm_hdo_ice))
724	call write_output("pdqvapD","pdqvapD apres microphysique"
725	& ,"kg/kg.s-1",pdqcloud(:,:,
726	& igcm_hdo_vap))
727	endif
728	call write_output("pdqccn2","pdqcloudccn apres microphysique"
729	& ,"kg/kg.s-1",pdqcloud(:,:,
730	& igcm_ccn_mass))
731	call write_output("pdqccnN2","pdqcloudccnN apres "//
732	& "microphysique","nb/kg.s-1",pdqcloud(:,:,
733	& igcm_ccn_number))
734	call write_output("pdqdust2", "pdqclouddust apres "//
735	& "microphysique","kg/kg.s-1",pdqcloud(:,:,
736	& igcm_dust_mass))
737	call write_output("pdqdustN2", "pdqclouddustN apres "//
738	& "microphysique","nb/kg.s-1",pdqcloud(:,:,
739	& igcm_dust_number))
740	c=======================================================================
741	#endif
742
743	END SUBROUTINE watercloud
744
745	subroutine ini_watercloud_mod(ngrid,nlayer,nq)
746	implicit none
747
748	integer,intent(in) :: ngrid ! number of atmospheric columns
749	integer,intent(in) :: nlayer ! number of atmospheric layers
750	integer,intent(in) :: nq ! number of tracers
751
752	allocate(zdqcloud(ngrid,nlayer,nq))
753	zdqcloud(:,:,:)=0
754	allocate(zdqscloud(ngrid,nq))
755	zdqscloud(:,:)=0
756
757	end subroutine ini_watercloud_mod
758
759
760	subroutine end_watercloud_mod
761	implicit none
762
763	if (allocated(zdqcloud)) deallocate(zdqcloud)
764	if (allocated(zdqscloud)) deallocate(zdqscloud)
765
766	end subroutine end_watercloud_mod
767
768	END MODULE watercloud_mod

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

Download in other formats:

Original Format