Context Navigation

watercloud_mod.F @ 3807

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

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

Download in other formats:

Original Format