Context Navigation

watercloud_mod.F @ 3118

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

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

Download in other formats:

Original Format