Context Navigation

watercloud_mod.F @ 3249

Last change on this file since 3249 was 2984, checked in by jnaar, 3 years 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