Context Navigation

watercloud.F @ 695

Last change on this file since 695 was 654, checked in by tnavarro, 13 years ago
bug correction if all dust disappears in cloud formation
File size: 15.4 KB

Rev	Line
[633]	1	SUBROUTINE watercloud(ngrid,nlay,ptimestep,
	2	& pplev,pplay,pdpsrf,pzlay,pt,pdt,
[626]	3	& pq,pdq,pdqcloud,pdtcloud,
[358]	4	& nq,tau,tauscaling,rdust,rice,nuice,
	5	& rsedcloud,rhocloud)
[633]	6	! to use 'getin'
	7	USE ioipsl_getincom
[38]	8	IMPLICIT NONE
	9
[633]	10
[38]	11	c=======================================================================
[358]	12	c Water-ice cloud formation
	13	c
	14	c Includes two different schemes:
	15	c - A simplified scheme (see simpleclouds.F)
	16	c - An improved microphysical scheme (see improvedclouds.F)
[38]	17	c
[633]	18	c There is a time loop specific to cloud formation
	19	c due to timescales smaller than the GCM integration timestep.
	20	c
[358]	21	c Authors: Franck Montmessin, Francois Forget, Ehouarn Millour,
[522]	22	c J.-B. Madeleine, Thomas Navarro
[38]	23	c
[633]	24	c 2004 - 2012
[38]	25	c=======================================================================
	26
	27	c-----------------------------------------------------------------------
	28	c declarations:
	29	c -------------
	30
	31	#include "dimensions.h"
	32	#include "dimphys.h"
	33	#include "comcstfi.h"
	34	#include "callkeys.h"
	35	#include "tracer.h"
	36	#include "comgeomfi.h"
[358]	37	#include "dimradmars.h"
[38]	38
	39	c Inputs:
	40	c ------
	41
	42	INTEGER ngrid,nlay
[633]	43	INTEGER nq ! nombre de traceurs
[38]	44	REAL ptimestep ! pas de temps physique (s)
	45	REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
	46	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
[633]	47	REAL pdpsrf(ngrid) ! tendence surf pressure
[38]	48	REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers
	49	REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K)
[633]	50	REAL pdt(ngrid,nlay) ! tendence temperature des autres param.
[38]	51
	52	real pq(ngrid,nlay,nq) ! traceur (kg/kg)
[633]	53	real pdq(ngrid,nlay,nq) ! tendence avant condensation (kg/kg.s-1)
[38]	54
[633]	55	REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
	56	REAL tauscaling(ngridmx) ! Convertion factor for dust amount
	57	real rdust(ngridmx,nlay) ! Dust geometric mean radius (m)
[38]	58
	59	c Outputs:
	60	c -------
	61
[633]	62	real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
	63	REAL pdtcloud(ngrid,nlay) ! tendence temperature due
	64	! a la chaleur latente
[38]	65
	66	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
	67	! (r_c in montmessin_2004)
	68	REAL nuice(ngrid,nlay) ! Estimated effective variance
	69	! of the size distribution
[633]	70	real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius
	71	real rhocloud(ngridmx,nlay) ! Cloud density (kg.m-3)
[38]	72
	73	c local:
	74	c ------
[633]	75
	76	! for ice radius computation
	77	REAL Mo,No
	78	REAl ccntyp
	79
	80	! for time loop
	81	INTEGER microstep ! time subsampling step variable
	82	INTEGER imicro ! time subsampling for coupled water microphysics & sedimentation
	83	SAVE imicro
	84	REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
	85	SAVE microtimestep
	86
	87	! tendency given by clouds (inside the micro loop)
	88	REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud
	89	REAL subpdtcloud(ngrid,nlay) ! cf. pdtcloud
[38]	90
[633]	91	! global tendency (clouds+physics)
	92	REAL subpdq(ngrid,nlay,nq) ! cf. pdqcloud
	93	REAL subpdt(ngrid,nlay) ! cf. pdtcloud
	94
	95	REAL CBRT
	96	EXTERNAL CBRT
	97
	98
	99	INTEGER iq,ig,l
[38]	100	LOGICAL,SAVE :: firstcall=.true.
	101
	102	c ** un petit test de coherence
	103	c --------------------------
	104
	105	IF (firstcall) THEN
	106	IF(ngrid.NE.ngridmx) THEN
	107	PRINT*,'STOP dans watercloud'
	108	PRINT*,'probleme de dimensions :'
	109	PRINT*,'ngrid =',ngrid
	110	PRINT*,'ngridmx =',ngridmx
	111	STOP
	112	ENDIF
	113
	114	if (nq.gt.nqmx) then
	115	write(,) 'stop in watercloud (nq.gt.nqmx)!'
	116	write(,) 'nq=',nq,' nqmx=',nqmx
	117	stop
	118	endif
	119
[358]	120	write(,) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
	121	write(,) " igcm_h2o_ice=",igcm_h2o_ice
[633]	122
	123	write(,) "time subsampling for microphysic ?"
	124	#ifdef MESOSCALE
	125	imicro = 2
	126	#else
	127	imicro = 15
	128	#endif
	129	call getin("imicro",imicro)
	130	write(,)"imicro = ",imicro
	131
	132	microtimestep = ptimestep/real(imicro)
	133	write(,)"Physical timestep is",ptimestep
	134	write(,)"Microphysics timestep is",microtimestep
[38]	135
	136	firstcall=.false.
	137	ENDIF ! of IF (firstcall)
[522]	138
[633]	139	c-----Initialization
	140	subpdq(1:ngrid,1:nlay,1:nq) = 0
	141	subpdt(1:ngrid,1:nlay) = 0
	142
	143	! default value if no ice
	144	rhocloud(1:ngrid,1:nlay) = rho_dust
[38]	145
[633]	146
	147
	148	c------------------------------------------------------------------
	149	c Time subsampling for microphysics
	150	c------------------------------------------------------------------
	151	DO microstep=1,imicro
[522]	152
[633]	153	c-------------------------------------------------------------------
	154	c 1. Tendencies:
	155	c------------------
[38]	156
[633]	157
	158	c------ Temperature tendency subpdt
	159	! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
	160	! If imicro=1 subpdt is the same as pdt
	161	DO l=1,nlay
	162	DO ig=1,ngrid
	163	subpdt(ig,l) = subpdt(ig,l)
	164	& + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
	165	ENDDO
	166	ENDDO
	167	c------ Tracers tendencies subpdq
	168	c------ At each micro timestep we add pdq in order to have a stepped entry
	169	IF (microphys) THEN
	170	DO l=1,nlay
	171	DO ig=1,ngrid
	172	subpdq(ig,l,igcm_dust_mass) =
	173	& subpdq(ig,l,igcm_dust_mass)
	174	& + pdq(ig,l,igcm_dust_mass)
	175	subpdq(ig,l,igcm_dust_number) =
	176	& subpdq(ig,l,igcm_dust_number)
	177	& + pdq(ig,l,igcm_dust_number)
	178	subpdq(ig,l,igcm_ccn_mass) =
	179	& subpdq(ig,l,igcm_ccn_mass)
	180	& + pdq(ig,l,igcm_ccn_mass)
	181	subpdq(ig,l,igcm_ccn_number) =
	182	& subpdq(ig,l,igcm_ccn_number)
	183	& + pdq(ig,l,igcm_ccn_number)
	184	ENDDO
	185	ENDDO
	186	ENDIF
	187	DO l=1,nlay
	188	DO ig=1,ngrid
	189	subpdq(ig,l,igcm_h2o_ice) =
	190	& subpdq(ig,l,igcm_h2o_ice)
	191	& + pdq(ig,l,igcm_h2o_ice)
	192	subpdq(ig,l,igcm_h2o_vap) =
	193	& subpdq(ig,l,igcm_h2o_vap)
	194	& + pdq(ig,l,igcm_h2o_vap)
	195	ENDDO
	196	ENDDO
	197
	198
	199	c-------------------------------------------------------------------
	200	c 2. Main call to the different cloud schemes:
	201	c------------------------------------------------
	202	IF (microphys) THEN
	203	CALL improvedclouds(ngrid,nlay,microtimestep,
	204	& pplay,pt,subpdt,
	205	& pq,subpdq,subpdqcloud,subpdtcloud,
	206	& nq,tauscaling)
	207
	208	ELSE
	209	CALL simpleclouds(ngrid,nlay,microtimestep,
	210	& pplay,pzlay,pt,subpdt,
	211	& pq,subpdq,subpdqcloud,subpdtcloud,
[645]	212	& nq,tau,rice)
[633]	213	ENDIF
	214
	215
	216	c-------------------------------------------------------------------
	217	c 3. Updating tendencies after cloud scheme:
	218	c-----------------------------------------------
	219
	220	IF (microphys) THEN
	221	DO l=1,nlay
	222	DO ig=1,ngrid
	223	subpdq(ig,l,igcm_dust_mass) =
	224	& subpdq(ig,l,igcm_dust_mass)
	225	& + subpdqcloud(ig,l,igcm_dust_mass)
	226	subpdq(ig,l,igcm_dust_number) =
	227	& subpdq(ig,l,igcm_dust_number)
	228	& + subpdqcloud(ig,l,igcm_dust_number)
	229	subpdq(ig,l,igcm_ccn_mass) =
	230	& subpdq(ig,l,igcm_ccn_mass)
	231	& + subpdqcloud(ig,l,igcm_ccn_mass)
	232	subpdq(ig,l,igcm_ccn_number) =
	233	& subpdq(ig,l,igcm_ccn_number)
	234	& + subpdqcloud(ig,l,igcm_ccn_number)
	235	ENDDO
	236	ENDDO
	237	ENDIF
	238	DO l=1,nlay
	239	DO ig=1,ngrid
	240	subpdq(ig,l,igcm_h2o_ice) =
	241	& subpdq(ig,l,igcm_h2o_ice)
	242	& + subpdqcloud(ig,l,igcm_h2o_ice)
	243	subpdq(ig,l,igcm_h2o_vap) =
	244	& subpdq(ig,l,igcm_h2o_vap)
	245	& + subpdqcloud(ig,l,igcm_h2o_vap)
	246	ENDDO
	247	ENDDO
	248
	249
	250	ENDDO ! of DO microstep=1,imicro
	251
	252	c-------------------------------------------------------------------
	253	c 6. Compute final tendencies after time loop:
	254	c------------------------------------------------
	255
	256	c------ Temperature tendency pdtcloud
	257	DO l=1,nlay
	258	DO ig=1,ngrid
	259	pdtcloud(ig,l) =
	260	& subpdt(ig,l)/real(imicro)-pdt(ig,l)
	261	ENDDO
	262	ENDDO
	263	c------ Tracers tendencies pdqcloud
	264	DO iq=1,nq
	265	DO l=1,nlay
	266	DO ig=1,ngrid
	267	pdqcloud(ig,l,iq) = subpdq(ig,l,iq)/real(imicro)
	268	& - pdq(ig,l,iq)
	269	ENDDO
	270	ENDDO
	271	ENDDO
	272
	273	c------- Due to stepped entry, other processes tendencies can add up to negative values
	274	c------- Therefore, enforce positive values and conserve mass
	275
	276	IF(microphys) THEN
	277	DO l=1,nlay
	278	DO ig=1,ngrid
[654]	279	IF ((pq(ig,l,igcm_dust_number) +
	280	& ptimestep* (pdq(ig,l,igcm_dust_number) +
	281	& pdqcloud(ig,l,igcm_dust_number)) .le. 1.)
	282	& .or. (pq(ig,l,igcm_dust_mass) +
	283	& ptimestep* (pdq(ig,l,igcm_dust_mass) +
	284	& pdqcloud(ig,l,igcm_dust_mass)) .le. 1.e-20)) THEN
	285	pdqcloud(ig,l,igcm_dust_number) =
	286	& - pq(ig,l,igcm_dust_number)/ptimestep
	287	& - pdq(ig,l,igcm_dust_number) + 1.
	288	pdqcloud(ig,l,igcm_ccn_number) =
	289	& -pdqcloud(ig,l,igcm_dust_number)
	290	pdqcloud(ig,l,igcm_dust_mass) =
	291	& - pq(ig,l,igcm_dust_mass)/ptimestep
	292	& - pdq(ig,l,igcm_dust_mass) + 1.e-20
	293	pdqcloud(ig,l,igcm_ccn_mass) =
	294	& -pdqcloud(ig,l,igcm_dust_mass)
	295	ENDIF
	296	IF ((pq(ig,l,igcm_ccn_number) +
[633]	297	& ptimestep* (pdq(ig,l,igcm_ccn_number) +
[654]	298	& pdqcloud(ig,l,igcm_ccn_number)) .le. 1.)
	299	& .or. (pq(ig,l,igcm_ccn_mass) +
	300	& ptimestep* (pdq(ig,l,igcm_ccn_mass) +
	301	& pdqcloud(ig,l,igcm_ccn_mass)) .le. 1.e-20)) THEN
[633]	302	pdqcloud(ig,l,igcm_ccn_number) =
	303	& - pq(ig,l,igcm_ccn_number)/ptimestep
[654]	304	& - pdq(ig,l,igcm_ccn_number) + 1.
[633]	305	pdqcloud(ig,l,igcm_dust_number) =
	306	& -pdqcloud(ig,l,igcm_ccn_number)
	307	pdqcloud(ig,l,igcm_ccn_mass) =
	308	& - pq(ig,l,igcm_ccn_mass)/ptimestep
[654]	309	& - pdq(ig,l,igcm_ccn_mass) + 1.e-20
[633]	310	pdqcloud(ig,l,igcm_dust_mass) =
	311	& -pdqcloud(ig,l,igcm_ccn_mass)
	312	ENDIF
	313	ENDDO
	314	ENDDO
	315	ENDIF
	316
	317	DO l=1,nlay
	318	DO ig=1,ngrid
	319	IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
	320	& (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice))
	321	& .le. 1.e-8) THEN
	322	pdqcloud(ig,l,igcm_h2o_ice) =
	323	& - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
	324	pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
	325	ENDIF
	326	IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
	327	& (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
	328	& .le. 1.e-8) THEN
	329	pdqcloud(ig,l,igcm_h2o_vap) =
	330	& - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
	331	pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
	332	ENDIF
	333	ENDDO
	334	ENDDO
	335
	336
	337	c------Update the ice and dust particle size "rice" for output or photochemistry
	338	c------Only rsedcloud is used for the water cycle
	339
	340	IF(scavenging) THEN
	341	DO l=1, nlay
	342	DO ig=1,ngrid
	343
	344	rdust(ig,l)=
	345	& CBRT(r3n_q*(pq(ig,l,igcm_dust_mass)+
	346	& (pdq(ig,l,igcm_dust_mass)
	347	& +pdqcloud(ig,l,igcm_dust_mass))*ptimestep)/
	348	& (pq(ig,l,igcm_dust_number)+
	349	& (pdq(ig,l,igcm_dust_number)+
	350	& pdqcloud(ig,l,igcm_dust_number)*ptimestep)))
	351	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
	352
	353	Mo = pq(ig,l,igcm_h2o_ice)
	354	& + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep
	355	& + (pq(ig,l,igcm_ccn_mass)
	356	& + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
	357	& *tauscaling(ig) + 1.e-30
	358	No = (pq(ig,l,igcm_ccn_number)
	359	& + pdqcloud(ig,l,igcm_ccn_number)*ptimestep)
	360	& *tauscaling(ig) + 1.e-30
	361	rhocloud(ig,l) = (pq(ig,l,igcm_h2o_ice) +
	362	& pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)
	363	& / Mo * rho_ice
	364	& + (pq(ig,l,igcm_ccn_mass)
	365	& + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
	366	& tauscaling(ig)/ Mo rho_dust
	367	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
	368	if ((Mo.lt.1.e-15) .or. (No.le.1.)) then
	369	rice(ig,l) = 1.e-8
	370	else
	371	!! AS: only perform computations if conditions not met
	372	rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.)
	373	rice(ig,l)=min(max(rice(ig,l),1.e-10),500.e-6)
	374	endif
	375	rsedcloud(ig,l)=max( rice(ig,l)(1.+nuice_sed)
	376	& (1.+nuice_sed)*(1.+nuice_sed),
	377	& rdust(ig,l) )
	378	rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
	379	ENDDO
	380	ENDDO
	381	ELSE
[645]	382
	383	if ((igcm_dust_mass.ne.0).and.(igcm_dust_number.ne.0)) then
	384	! recompute rdust(), if we have dust_mass & dust_number tracers
	385	DO l=1,nlay
[633]	386	DO ig=1,ngrid
[645]	387	rdust(ig,l)=
[633]	388	& CBRT(r3n_q*(pq(ig,l,igcm_dust_mass)+
	389	& (pdq(ig,l,igcm_dust_mass)
	390	& +pdqcloud(ig,l,igcm_dust_mass))*ptimestep)/
	391	& (pq(ig,l,igcm_dust_number)+
	392	& (pdq(ig,l,igcm_dust_number)+
	393	& pdqcloud(ig,l,igcm_dust_number)*ptimestep)))
[645]	394	rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
	395	ENDDO
	396	ENDDO
	397	endif ! of if ((igcm_dust_mass.ne.0).and.(igcm_dust_number.ne.0))
	398
	399	DO l=1,nlay
	400	DO ig=1,ngrid
[633]	401	ccntyp =
	402	& 1.3e+8max(tau(ig,1),0.001)/0.1exp(-pzlay(ig,l)/10000.)
	403	ccntyp = ccntyp /ccn_factor
	404	rice(ig,l)=max( CBRT ( ((pq(ig,l,igcm_h2o_ice)
	405	& + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)/rho_ice
	406	& +ccntyp(4./3.)pirdust(ig,l)rdust(ig,l)*rdust(ig,l))
	407	& /(ccntyp4./3.pi) ), rdust(ig,l))
	408	rsedcloud(ig,l)=max( rice(ig,l)(1.+nuice_sed)
	409	& (1.+nuice_sed)*(1.+nuice_sed),
	410	& rdust(ig,l) )
	411	rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
	412	ENDDO
	413	ENDDO
	414	ENDIF ! of IF(scavenging)
	415
	416
	417	! used for rad. transfer calculations
	418	! nuice is constant because a lognormal distribution is prescribed
	419	nuice(1:ngrid,1:nlay)=nuice_ref
	420
	421
	422	!--------------------------------------------------------------
	423	!--------------------------------------------------------------
	424
	425
[358]	426	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
	427	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	428	c Then that should not affect the ice particle radius
	429	do ig=1,ngridmx
	430	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
	431	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
	432	& rice(ig,2)=rice(ig,3)
	433	rice(ig,1)=rice(ig,2)
	434	end if
	435	end do
[38]	436
[633]	437	c=======================================================================
	438
[38]	439	END
	440

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

Download in other formats:

Original Format