Context Navigation

watercloud_mod.F @ 2937

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

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

Download in other formats:

Original Format