Context Navigation

watercloud_mod.F @ 2966

Last change on this file since 2966 was 2966, checked in by jnaar, 19 months ago
Architecture changes in watercloud_mod, improvedclouds_mod in preparation for adaptative timestep of clouds (not yet implemented). Changes prevent bit by bit correspondance with previous revision. "simpleclouds" routine broken.
File size: 31.0 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, adding different subtimestep on each grid cell
50	c plus not doing microphysics if no water present
51	c plus simpleclouds no longer in the loop for microphysics
52	c=======================================================================
53
54	c-----------------------------------------------------------------------
55	c declarations:
56	c -------------
57
58	include "callkeys.h"
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 :: computed_micro(ngrid,nlay) ! Check if microphy was done inthis cell (logical)
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
158
159	c ** un petit test de coherence
160	c --------------------------
161
162	IF (firstcall) THEN
163
164	if (nq.gt.nqmx) then
165	write(,) 'stop in watercloud (nq.gt.nqmx)!'
166	write(,) 'nq=',nq,' nqmx=',nqmx
167	call abort_physic("watercloud","nq.gt.nqmx",1)
168	endif
169
170	write(,) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
171	write(,) " igcm_h2o_ice=",igcm_h2o_ice
172
173	write(,) "time subsampling for microphysic ?"
174	#ifdef MESOSCALE
175	imicro = 2
176	#else
177	imicro = 30
178	#endif
179	call getin_p("imicro",imicro)
180	write(,)"watercloud: imicro = ",imicro
181
182	firstcall=.false.
183	ENDIF ! of IF (firstcall)
184
185	!! AS: moved out of firstcall to allow nesting+evoluting timestep
186	!! TBD: consider possible diff imicro with domains?
187	microtimestep = ptimestep/real(imicro)
188	if (microtimestep/=microtimestep_prev) then
189	! only tell the world if microtimestep has changed
190	write(,)"watercloud: Physical timestep is ",ptimestep
191	write(,)"watercloud: Microphysics timestep is ",microtimestep
192	microtimestep_prev=microtimestep
193	endif
194
195	c-----Initialization
196	sum_subpdq(1:ngrid,1:nlay,1:nq) = 0
197	sum_subpdt(1:ngrid,1:nlay) = 0
198
199	! default value if no ice
200	rhocloud(1:ngrid,1:nlay) = rho_dust
201
202	c-------------------------------------------------------------------
203	c 0. Representation of sub-grid water ice clouds
204	c------------------
205	c-----Initialization
206	pteff(1:ngrid,1:nlay) = 0
207	pqeff(1:ngrid,1:nlay,1:nq) = 0
208	DO l=1,nlay
209	DO ig=1,ngrid
210	pteff(ig,l)=pt(ig,l)
211	END DO
212	END DO
213	DO l=1,nlay
214	DO ig=1,ngrid
215	DO iq=1,nq
216	pqeff(ig,l,iq)=pq(ig,l,iq)
217	ENDDO
218	ENDDO
219	ENDDO
220	c-----Tendencies
221	DO l=1,nlay
222	DO ig=1,ngrid
223	ztclf(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
224	ENDDO
225	ENDDO
226	DO l=1,nlay
227	DO ig=1,ngrid
228	DO iq=1,nq
229	zqclf(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep
230	ENDDO
231	ENDDO
232	ENDDO
233	c-----Effective temperature calculation
234	IF (CLFvarying) THEN
235	spant=3.0 ! delta T for the temprature distribution
236	mincloud=0.1 ! min cloudfrac when there is ice
237	IF (flagcloud) THEN
238	write(,) "Delta T", spant
239	write(,) "mincloud", mincloud
240	flagcloud=.false.
241	END IF
242	!CALL watersat(ngrid*nlay,ztclf,pplay,zqsat) !MV17: we dont need zqsat in the CLFvarying scheme
243	zqvap=zqclf(:,:,igcm_h2o_vap)
244	zqice=zqclf(:,:,igcm_h2o_ice)
245	CALL tcondwater(ngrid*nlay,pplay,zqvap+zqice,tcond)
246	DO l=1,nlay
247	DO ig=1,ngrid
248	zdelt=spant !MAX(spant*ztclf(ig,l),1.e-12), now totally in K. Fixed width
249	IF (tcond(ig,l) .ge. (ztclf(ig,l)+zdelt)) THEN
250	pteff(ig,l)=ztclf(ig,l)
251	cloudfrac(ig,l)=1.
252	ELSE IF (tcond(ig,l) .le. (ztclf(ig,l)-zdelt)) THEN
253	pteff(ig,l)=ztclf(ig,l)-zdelt
254	cloudfrac(ig,l)=mincloud
255	ELSE
256	cloudfrac(ig,l)=(tcond(ig,l)-ztclf(ig,l)+zdelt)/
257	& (2.0*zdelt)
258	pteff(ig,l)=(tcond(ig,l)+ztclf(ig,l)-zdelt)/2.
259	END IF
260	pteff(ig,l)=pteff(ig,l)-pdt(ig,l)*ptimestep
261	IF (cloudfrac(ig,l).le.mincloud) THEN !MV17: replaced .le.0 by .le.mincloud
262	cloudfrac(ig,l)=mincloud
263	ELSE IF (cloudfrac(ig,l).gt.1) THEN
264	cloudfrac(ig,l)=1.
265	END IF
266	ENDDO
267	ENDDO
268	c-----Calculation of the total cloud coverage of the column
269	DO ig=1,ngrid
270	totcloudfrac(ig) = 0.
271	norm=0.
272	DO l=1,nlay
273	ponder=zqice(ig,l)*(pplev(ig,l) - pplev(ig,l+1))
274	totcloudfrac(ig) = totcloudfrac(ig)
275	& + cloudfrac(ig,l)*ponder
276	norm=norm+ponder
277	ENDDO
278	totcloudfrac(ig)=MAX(totcloudfrac(ig)/norm,1.e-12) ! min value if NaNs
279	ENDDO
280	c-----Effective tracers quantities in the cloud fraction
281	IF (microphys) THEN
282	pqeff(:,:,igcm_ccn_mass)=pq(:,:,igcm_ccn_mass)/
283	& cloudfrac(:,:)
284	pqeff(:,:,igcm_ccn_number)=pq(:,:,igcm_ccn_number)/
285	& cloudfrac(:,:)
286	END IF ! end if (microphys)
287	pqeff(:,:,igcm_h2o_ice)=pq(:,:,igcm_h2o_ice)/
288	& cloudfrac(:,:)
289	!! CLFvarying outputs
290	! CALL write_output('pqeffice','pqeffice',
291	! & 'kg/kg',pqeff(:,:,igcm_h2o_ice))
292	! CALL write_output('pteff','pteff',
293	! & 'K',pteff(:,:))
294	! CALL write_output('tcond','tcond',
295	! & 'K',tcond(:,:))
296	! CALL write_output('cloudfrac','cloudfrac',
297	! & 'K',cloudfrac(:,:))
298	END IF ! end if (CLFvarying)
299	c------------------------------------------------------------------
300	c Time subsampling for microphysics
301	c------------------------------------------------------------------
302	rhocloud(1:ngrid,1:nlay) = rho_dust
303
304
305	c Initialisation of all the stuff JN2023
306	c computed_micro(1:ngrid,1:nlay)=.false.
307	computed_micro(1:ngrid,1:nlay)=0.
308	zt_micro(:,:)=pt(:,:)
309	zq_micro(:,:,:)=pq(:,:,:)
310	zq_micro(:,:,:)=pq(:,:,:)
311	call watersat(ngrid*nlay,zt_micro,pplay,zqsat_micro)
312	zpotcond_inst=zq_micro(:,:,igcm_h2o_vap) - zqsat_micro
313	call watersat(ngridnlay,zt_micro+pdtptimestep,pplay,zqsat_micro)
314	zpotcond_full=(zq_micro(:,:,igcm_h2o_vap)+
315	& pdq(:,:,igcm_h2o_vap)*ptimestep) - zqsat_micro
316	zimicro(1:ngrid,1:nlay)=imicro
317	if (cloud_adapt_ts) then
318	call adapt_imicro(ptimestep,zpotcond(ig,l),
319	$ zimicro(ig,l))
320	endif! (cloud_adapt_ts) then
321	DO l=1,nlay
322	DO ig=1,ngrid
323	c Start by computing the condensable water vapor amount
324	if (zpotcond_full(ig,l).gt.0.) then
325	zpotcond(ig,l)=max(zpotcond_inst(ig,l),zpotcond_full(ig,l))
326	else if (zpotcond_full(ig,l).le.0.) then
327	zpotcond(ig,l)=min(zpotcond_inst(ig,l),zpotcond_full(ig,l))
328	endif! (zpotcond_full.gt.0.) then
329	microtimestep=ptimestep/real(zimicro(ig,l))
330	c Check if microphysics is even needed, that is if enough action
331	c is happening water-wise
332	if ((pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
333	& .gt.1e-22) .or. (abs(zpotcond(ig,l)).gt.1e-22)) then
334	c Eventuellement sortir simpleclouds de la boucle egalement
335	computed_micro(ig,l)=1.
336	DO microstep=1,zimicro(ig,l)
337
338
339	c JN : incrementing after main microphysics scheme
340	c Previously we were incrementing tendencies, we now
341	c increment tracers and temperature directly
342	c We are thus starting at the end of the first iteration
343	c
344	c-------------------------------------------------------------------
345	c 1. Main call to the different cloud schemes:
346	c------------------------------------------------
347	IF (microphys) THEN
348	CALL improvedclouds(microtimestep,
349	& pplay(ig,l),zt_micro(ig,l),
350	& zq_micro(ig,l,:),subpdqcloud(ig,l,:),
351	& subpdtcloud(ig,l),nq,tauscaling(ig),mmean(ig,l))
352
353	ELSE
354	c Simpleclouds should maybe be taken out and put in a specific loop ?
355	CALL simpleclouds(ngrid,nlay,microtimestep,
356	& pplay,pzlay,pteff,sum_subpdt,
357	& pqeff,sum_subpdq,subpdqcloud,subpdtcloud,
358	& nq,tau,rice)
359	ENDIF
360
361	c-------------------------------------------------------------------
362	c 2. Updating tracers and temperature after cloud scheme:
363	c-----------------------------------------------
364
365	IF (microphys) THEN
366	zq_micro(ig,l,igcm_dust_mass) =
367	& zq_micro(ig,l,igcm_dust_mass)+(pdq(ig,l,igcm_dust_mass)
368	& +subpdqcloud(ig,l,igcm_dust_mass))*microtimestep
369	zq_micro(ig,l,igcm_dust_number) =
370	& zq_micro(ig,l,igcm_dust_number)
371	& +(pdq(ig,l,igcm_dust_number)
372	& + subpdqcloud(ig,l,igcm_dust_number))*microtimestep
373	zq_micro(ig,l,igcm_ccn_mass) =
374	& zq_micro(ig,l,igcm_ccn_mass) +
375	& (pdq(ig,l,igcm_ccn_mass)
376	& +subpdqcloud(ig,l,igcm_ccn_mass))*microtimestep
377	zq_micro(ig,l,igcm_ccn_number) =
378	& zq_micro(ig,l,igcm_ccn_number) +
379	& (pdq(ig,l,igcm_ccn_number)
380	& + subpdqcloud(ig,l,igcm_ccn_number))*microtimestep
381	ENDIF
382	zq_micro(ig,l,igcm_h2o_ice) =
383	& zq_micro(ig,l,igcm_h2o_ice)+
384	& (pdq(ig,l,igcm_h2o_ice)
385	& + subpdqcloud(ig,l,igcm_h2o_ice))*microtimestep
386	zq_micro(ig,l,igcm_h2o_vap) =
387	& zq_micro(ig,l,igcm_h2o_vap)+
388	& (pdq(ig,l,igcm_h2o_vap)
389	& + subpdqcloud(ig,l,igcm_h2o_vap))*microtimestep
390
391	IF (hdo) THEN
392	zq_micro(ig,l,igcm_hdo_ice) =
393	& zq_micro(ig,l,igcm_hdo_ice)+
394	& (pdq(ig,l,igcm_hdo_ice)
395	& + subpdqcloud(ig,l,igcm_hdo_ice))*microtimestep
396	zq_micro(ig,l,igcm_hdo_vap) =
397	& zq_micro(ig,l,igcm_hdo_vap)+
398	& (pdq(ig,l,igcm_hdo_vap)
399	& + subpdqcloud(ig,l,igcm_hdo_vap))*microtimestep
400	ENDIF ! hdo
401
402	c Could also set subpdtcloud to 0 if not activice to make it simpler
403	zt_micro(ig,l) = zt_micro(ig,l)+
404	& pdt(ig,l)*microtimestep
405	IF (activice) THEN
406	zt_micro(ig,l) = zt_micro(ig,l)+
407	& subpdtcloud(ig,l)*microtimestep
408	ENDIF
409	! !! Example of how to use writediagmicrofi useful to
410	! !! get outputs at each microphysical sub-timestep (better to be used in 1D)
411	! CALL WRITEDIAGMICROFI(ngrid,imicro,microstep,
412	! & microtimestep,'subpdtcloud',
413	! & 'subpdtcloud','K/s',1,subpdtcloud(:,:))
414
415	ENDDO ! of DO microstep=1,imicro
416	endif! (zq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
417	! & .gt.1e-22).or.(abs(zpotcond).gt.1e-22) then
418	ENDDO ! ig=1,ngrid
419	ENDDO ! l=1,nlay
420
421
422	c------ Useful outputs to check how it went
423	call write_output("computed_micro","computed_micro "//
424	& "after microphysics","logical",computed_micro(:,:))
425	call write_output("zimicro","Used number of subtimestep "//
426	& "in cloud microphysics","integer",real(zimicro(:,:)))
427	c-------------------------------------------------------------------
428	c 3. Compute final tendencies after time loop:
429	c------------------------------------------------
430
431	c------ Temperature tendency pdtcloud
432	DO l=1,nlay
433	DO ig=1,ngrid
434	pdtcloud(ig,l) = -pdt(ig,l)+
435	& (zt_micro(ig,l)-pt(ig,l)) / ptimestep
436	ENDDO
437	ENDDO
438
439	c------ Tracers tendencies pdqcloud
440	DO l=1,nlay
441	DO ig=1,ngrid
442	pdqcloud(ig,l,igcm_h2o_ice) =
443	& -pdq(ig,l,igcm_h2o_ice)+
444	& (zq_micro(ig,l,igcm_h2o_ice) -
445	& pq(ig,l,igcm_h2o_ice)) / ptimestep
446	pdqcloud(ig,l,igcm_h2o_vap) =
447	& -pdq(ig,l,igcm_h2o_vap)+
448	& (zq_micro(ig,l,igcm_h2o_vap) -
449	& pq(ig,l,igcm_h2o_vap)) / ptimestep
450	IF (hdo) THEN
451	pdqcloud(ig,l,igcm_hdo_ice) =
452	& -pdq(ig,l,igcm_hdo_ice)+
453	& (zq_micro(ig,l,igcm_hdo_ice) -
454	& pq(ig,l,igcm_hdo_ice)) / ptimestep
455	pdqcloud(ig,l,igcm_hdo_vap) =
456	& -pdq(ig,l,igcm_hdo_vap)+
457	& (zq_micro(ig,l,igcm_hdo_vap) -
458	& pq(ig,l,igcm_hdo_vap)) / ptimestep
459	ENDIF !hdo
460	ENDDO
461	ENDDO
462
463	IF(microphys) THEN
464	DO l=1,nlay
465	DO ig=1,ngrid
466	pdqcloud(ig,l,igcm_ccn_mass) =
467	& -pdq(ig,l,igcm_ccn_mass)+
468	& (zq_micro(ig,l,igcm_ccn_mass) -
469	& pq(ig,l,igcm_ccn_mass)) / ptimestep
470	pdqcloud(ig,l,igcm_ccn_number) =
471	& -pdq(ig,l,igcm_ccn_number)+
472	& (zq_micro(ig,l,igcm_ccn_number) -
473	& pq(ig,l,igcm_ccn_number)) / ptimestep
474	ENDDO
475	ENDDO
476	ENDIF
477
478	IF(scavenging) THEN
479	DO l=1,nlay
480	DO ig=1,ngrid
481	pdqcloud(ig,l,igcm_dust_mass) =
482	& -pdq(ig,l,igcm_dust_mass)+
483	& (zq_micro(ig,l,igcm_dust_mass) -
484	& pq(ig,l,igcm_dust_mass)) / ptimestep
485	pdqcloud(ig,l,igcm_dust_number) =
486	& -pdq(ig,l,igcm_dust_number)+
487	& (zq_micro(ig,l,igcm_dust_number) -
488	& pq(ig,l,igcm_dust_number)) / ptimestep
489	ENDDO
490	ENDDO
491	ENDIF
492
493	c------- Due to stepped entry, other processes tendencies can add up to negative values
494	c------- Therefore, enforce positive values and conserve mass
495	IF(microphys) THEN
496	DO l=1,nlay
497	DO ig=1,ngrid
498	IF ((pq(ig,l,igcm_ccn_number) +
499	& ptimestep* (pdq(ig,l,igcm_ccn_number) +
500	& pdqcloud(ig,l,igcm_ccn_number)) .le. 1.)
501	& .or. (pq(ig,l,igcm_ccn_mass) +
502	& ptimestep* (pdq(ig,l,igcm_ccn_mass) +
503	& pdqcloud(ig,l,igcm_ccn_mass)) .le. 1.e-20)) THEN
504	pdqcloud(ig,l,igcm_ccn_number) =
505	& - pq(ig,l,igcm_ccn_number)/ptimestep
506	& - pdq(ig,l,igcm_ccn_number) + 1.
507	pdqcloud(ig,l,igcm_dust_number) =
508	& -pdqcloud(ig,l,igcm_ccn_number)
509	pdqcloud(ig,l,igcm_ccn_mass) =
510	& - pq(ig,l,igcm_ccn_mass)/ptimestep
511	& - pdq(ig,l,igcm_ccn_mass) + 1.e-20
512	pdqcloud(ig,l,igcm_dust_mass) =
513	& -pdqcloud(ig,l,igcm_ccn_mass)
514	ENDIF
515	ENDDO
516	ENDDO
517	ENDIF
518
519	IF(scavenging) THEN
520	DO l=1,nlay
521	DO ig=1,ngrid
522	IF ((pq(ig,l,igcm_dust_number) +
523	& ptimestep* (pdq(ig,l,igcm_dust_number) +
524	& pdqcloud(ig,l,igcm_dust_number)) .le. 1.)
525	& .or. (pq(ig,l,igcm_dust_mass) +
526	& ptimestep* (pdq(ig,l,igcm_dust_mass) +
527	& pdqcloud(ig,l,igcm_dust_mass)) .le. 1.e-20)) THEN
528	pdqcloud(ig,l,igcm_dust_number) =
529	& - pq(ig,l,igcm_dust_number)/ptimestep
530	& - pdq(ig,l,igcm_dust_number) + 1.
531	pdqcloud(ig,l,igcm_ccn_number) =
532	& -pdqcloud(ig,l,igcm_dust_number)
533	pdqcloud(ig,l,igcm_dust_mass) =
534	& - pq(ig,l,igcm_dust_mass)/ptimestep
535	& - pdq(ig,l,igcm_dust_mass) + 1.e-20
536	pdqcloud(ig,l,igcm_ccn_mass) =
537	& -pdqcloud(ig,l,igcm_dust_mass)
538	ENDIF
539	ENDDO
540	ENDDO
541	ENDIF
542
543	DO l=1,nlay
544	DO ig=1,ngrid
545
546	IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
547	& (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice))
548	& .le. qperemin) THEN
549	pdqcloud(ig,l,igcm_h2o_ice) =
550	& - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
551	pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
552	if (hdo) then
553	pdqcloud(ig,l,igcm_hdo_ice) =
554	& - pq(ig,l,igcm_hdo_ice)/ptimestep - pdq(ig,l,igcm_hdo_ice)
555	pdqcloud(ig,l,igcm_hdo_vap) = -pdqcloud(ig,l,igcm_hdo_ice)
556	endif
557	ENDIF
558	IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
559	& (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
560	& .le. qperemin) THEN
561	pdqcloud(ig,l,igcm_h2o_vap) =
562	& - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
563	pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
564	if (hdo) then
565	pdqcloud(ig,l,igcm_hdo_vap) =
566	& - pq(ig,l,igcm_hdo_vap)/ptimestep - pdq(ig,l,igcm_hdo_vap)
567	pdqcloud(ig,l,igcm_hdo_ice) = -pdqcloud(ig,l,igcm_hdo_vap)
568	endif
569	ENDIF
570
571	ENDDO
572	ENDDO
573
574	c------Update the ice and dust particle size "rice" for output or photochemistry
575	c------Only rsedcloud is used for the water cycle
576
577	IF(scavenging) THEN
578	DO l=1, nlay
579	DO ig=1,ngrid
580
581	call updaterdust(
582	& pq(ig,l,igcm_dust_mass) + ! dust mass
583	& (pdq(ig,l,igcm_dust_mass) + ! dust mass
584	& pdqcloud(ig,l,igcm_dust_mass))*ptimestep, ! dust mass
585	& pq(ig,l,igcm_dust_number) + ! dust number
586	& (pdq(ig,l,igcm_dust_number) + ! dust number
587	& pdqcloud(ig,l,igcm_dust_number))*ptimestep, ! dust number
588	& rdust(ig,l))
589
590	ENDDO
591	ENDDO
592	ENDIF
593
594	IF(microphys) THEN
595
596	! In case one does not want to allow supersatured water when using microphysics.
597	! Not done by default.
598	IF(.not.supersat) THEN
599	! !! initial mixing ratios for initial D/H ratio calculation
600	zq0(:,:,:) = pq(:,:,:) + pdq(:,:,:)*ptimestep
601	zt = pt + (pdt+pdtcloud)*ptimestep
602	call watersat(ngrid*nlay,zt,pplay,zqsat)
603	DO l=1, nlay
604	DO ig=1,ngrid
605	IF (pq(ig,l,igcm_h2o_vap)
606	& + (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
607	& * ptimestep .ge. zqsat(ig,l)) THEN
608	pdqcloud(ig,l,igcm_h2o_vap) =
609	& (zqsat(ig,l) - pq(ig,l,igcm_h2o_vap))/ptimestep
610	& - pdq(ig,l,igcm_h2o_vap)
611	pdqcloud(ig,l,igcm_h2o_ice) =
612	& -pdqcloud(ig,l,igcm_h2o_vap)
613	! !! HDO ice flux has to be modified in consequence
614	IF (hdo) THEN
615	! !! Logically only condensation can happen in this case
616	IF (pdqcloud(ig,l,igcm_h2o_ice) .gt. 0.0) THEN
617	IF ( zq0(ig,l,igcm_h2o_vap) .gt. qperemin ) THEN
618	! !! Lamb et al. 2017
619	alpha(ig,l) = exp(13525./zt(ig,l)**2.-5.59d-2)
620	pdqcloud(ig,l,igcm_hdo_ice) =
621	& pdqcloud(ig,l,igcm_h2o_ice)alpha(ig,l)
622	& ( zq0(ig,l,igcm_hdo_vap)
623	& /zq0(ig,l,igcm_h2o_vap) )
624	pdqcloud(ig,l,igcm_hdo_ice) =
625	& min(pdqcloud(ig,l,igcm_hdo_ice),
626	& zq0(ig,l,igcm_hdo_vap)/ptimestep)
627	ELSE
628	pdqcloud(ig,l,igcm_hdo_ice) = 0.
629	ENDIF
630	!! sublimation
631	ELSE
632	IF ( zq0(ig,l,igcm_h2o_ice) .gt. qperemin ) THEN
633	pdqcloud(ig,l,igcm_hdo_ice) =
634	& pdqcloud(ig,l,igcm_h2o_ice)*
635	& ( zq0(ig,l,igcm_hdo_ice)
636	& /zq0(ig,l,igcm_h2o_ice) )
637	pdqcloud(ig,l,igcm_hdo_ice) =
638	& max(pdqcloud(ig,l,igcm_hdo_ice),
639	& -zq0(ig,l,igcm_hdo_ice)/ptimestep)
640	ELSE
641	pdqcloud(ig,l,igcm_hdo_ice) = 0.
642	ENDIF
643	ENDIF !IF (pdqcloud(ig,l,igcm_h2o_ice).gt.0.)
644	pdqcloud(ig,l,igcm_hdo_vap) =
645	& -pdqcloud(ig,l,igcm_hdo_ice)
646	ENDIF !IF (hdo)
647	! no need to correct ccn_number, updaterad can handle this properly.
648	ENDIF
649	ENDDO
650	ENDDO
651	ENDIF
652
653	DO l=1, nlay
654	DO ig=1,ngrid
655
656	call updaterice_micro(
657	& pq(ig,l,igcm_h2o_ice) + ! ice mass
658	& (pdq(ig,l,igcm_h2o_ice) + ! ice mass
659	& pdqcloud(ig,l,igcm_h2o_ice))*ptimestep, ! ice mass
660	& pq(ig,l,igcm_ccn_mass) + ! ccn mass
661	& (pdq(ig,l,igcm_ccn_mass) + ! ccn mass
662	& pdqcloud(ig,l,igcm_ccn_mass))*ptimestep, ! ccn mass
663	& pq(ig,l,igcm_ccn_number) + ! ccn number
664	& (pdq(ig,l,igcm_ccn_number) + ! ccn number
665	& pdqcloud(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
666	& tauscaling(ig),rice(ig,l),rhocloud(ig,l))
667
668	ENDDO
669	ENDDO
670
671	ELSE ! no microphys
672
673	DO l=1,nlay
674	DO ig=1,ngrid
675
676	call updaterice_typ(
677	& pq(ig,l,igcm_h2o_ice) + ! ice mass
678	& (pdq(ig,l,igcm_h2o_ice) + ! ice mass
679	& pdqcloud(ig,l,igcm_h2o_ice))*ptimestep, ! ice mass
680	& tau(ig,1),pzlay(ig,l),rice(ig,l))
681
682	ENDDO
683	ENDDO
684
685	ENDIF ! of IF(microphys)
686
687
688
689	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
690	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
691	c Then that should not affect the ice particle radius
692	do ig=1,ngrid
693	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
694	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
695	& rice(ig,2)=rice(ig,3)
696	rice(ig,1)=rice(ig,2)
697	end if
698	end do
699
700
701	DO l=1,nlay
702	DO ig=1,ngrid
703	rsedcloud(ig,l)=max(rice(ig,l)*
704	& (1.+nuice_sed)(1.+nuice_sed)(1.+nuice_sed),
705	& rdust(ig,l))
706	! rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
707	ENDDO
708	ENDDO
709
710	! used for rad. transfer calculations
711	! nuice is constant because a lognormal distribution is prescribed
712	nuice(1:ngrid,1:nlay)=nuice_ref
713
714	c------Update tendencies for sub-grid water ice clouds
715	IF (CLFvarying) THEN
716	DO ig=1,ngrid
717	DO l=1,nlay
718	pdqcloud(ig,l,igcm_dust_mass)=pdqcloud(ig,l,igcm_dust_mass)
719	& *cloudfrac(ig,l)
720	pdqcloud(ig,l,igcm_ccn_mass)=pdqcloud(ig,l,igcm_ccn_mass)
721	& *cloudfrac(ig,l)
722	pdqcloud(ig,l,igcm_dust_number)=pdqcloud(ig,l,
723	& igcm_dust_number) *cloudfrac(ig,l)
724	pdqcloud(ig,l,igcm_ccn_number)=pdqcloud(ig,l,
725	& igcm_ccn_number) *cloudfrac(ig,l)
726	pdqcloud(ig,l,igcm_h2o_vap)=pdqcloud(ig,l,
727	& igcm_h2o_vap) *cloudfrac(ig,l)
728	pdqcloud(ig,l,igcm_h2o_ice)=pdqcloud(ig,l,
729	& igcm_h2o_ice) *cloudfrac(ig,l)
730	ENDDO
731	ENDDO
732	pdtcloud(:,:)=pdtcloud(:,:)*cloudfrac(:,:)
733	ENDIF
734	#ifndef MESOSCALE
735	c=======================================================================
736	call write_output("watercloud_pdqh2oice","pdqcloud_h2o_ice "//
737	& "after microphysics","kg/kg.s-1",pdqcloud(:,:,igcm_h2o_ice))
738	call write_output("watercloud_pdqh2ovap","pdqcloud_h2o_vap "//
739	& "after microphysics","kg/kg.s-1",pdqcloud(:,:,igcm_h2o_vap))
740	if (hdo) then
741	call write_output("watercloud_pdqhdoice","pdqcloud_hdo_ice "//
742	& "after microphysics","kg/kg.s-1",pdqcloud(:,:,igcm_hdo_ice))
743	call write_output("watercloud_pdqhdovap","pdqcloud_hdo_vap "//
744	& "after microphysics","kg/kg.s-1",pdqcloud(:,:,igcm_hdo_vap))
745	endif
746	! call write_output("pdqccn2","pdqcloudccn apres microphysique"
747	! & ,"kg/kg.s-1",pdqcloud(:,:,
748	! & igcm_ccn_mass))
749	! call write_output("pdqccnN2","pdqcloudccnN apres "//
750	! & "microphysique","nb/kg.s-1",pdqcloud(:,:,
751	! & igcm_ccn_number))
752	! call write_output("pdqdust2", "pdqclouddust apres "//
753	! & "microphysique","kg/kg.s-1",pdqcloud(:,:,
754	! & igcm_dust_mass))
755	! call write_output("pdqdustN2", "pdqclouddustN apres "//
756	! & "microphysique","nb/kg.s-1",pdqcloud(:,:,
757	! & igcm_dust_number))
758	c=======================================================================
759	#endif
760
761	END SUBROUTINE watercloud
762
763	subroutine ini_watercloud_mod(ngrid,nlayer,nq)
764	implicit none
765
766	integer,intent(in) :: ngrid ! number of atmospheric columns
767	integer,intent(in) :: nlayer ! number of atmospheric layers
768	integer,intent(in) :: nq ! number of tracers
769
770	allocate(zdqcloud(ngrid,nlayer,nq))
771	zdqcloud(:,:,:)=0
772	allocate(zdqscloud(ngrid,nq))
773	zdqscloud(:,:)=0
774
775	end subroutine ini_watercloud_mod
776
777
778	subroutine end_watercloud_mod
779	implicit none
780
781	if (allocated(zdqcloud)) deallocate(zdqcloud)
782	if (allocated(zdqscloud)) deallocate(zdqscloud)
783
784	end subroutine end_watercloud_mod
785
786
787
788	SUBROUTINE adapt_imicro(ptimestep,potcond,
789	$ zimicro)
790
791	c Pas de temps adaptatif pour les nuages
792
793	real,intent(in) :: ptimestep ! total duration of physics (sec)
794	real,intent(in) :: potcond ! total duration of physics (sec)
795	real :: alpha, beta ! total duration of physics (sec)
796	integer,intent(out) :: zimicro ! number of ptimestep division
797
798	c zimicro = ptimestepalphapotcond**beta
799	zimicro = 30
800
801	END SUBROUTINE adapt_imicro
802
803
804	END MODULE watercloud_mod

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

Download in other formats:

Original Format