Context Navigation

watercloud_mod.F @ 2814

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

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

Download in other formats:

Original Format