Context Navigation

watercloud_mod.F @ 2541

Last change on this file since 2541 was 2516, checked in by mvals, 4 years ago

Mars GCM:
Correction of the threshold value used to prevent from negative values in watercloud_mod.F: the threshold is set to 1e-16 (=qperemin, an already
existing parameter for the water cycle), instead of 1e-8. This former value was at the origin of a strong and unrealistic drying out of the upper
atmosphere, when the photochemistry was active, during the first part of the martian year (months 3 and 4). (bug reported by Francisco
Gonzalez-Galindo).
MV

File size: 26.7 KB

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

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

Download in other formats:

Original Format