Context Navigation

watercloud_mod.F @ 2347

Last change on this file since 2347 was 2312, checked in by lrossi, 5 years ago

MARS GCM:
Implementation of HDO in the physics
New tracers hdo_vap and hdo_ice are added. Cf. traceur.def.isotopes in deftank for exemple of traceur.def.
All HDO related computations are activated by flag hdo=.true. in callphys.def. (see callphys.def.hdo in deftank for an example)
Additional option hdofrac (true by default) allows for turning on/off fractionation (for tests).
For now, only functional with simplified cloud scheme (so microphys=.false. and activice=.false. also recommended).
Initialisation of start files can be done with option 'inihdo' in newstart.
LR

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

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

Download in other formats:

Original Format