Context Navigation

watercloud_mod.F @ 2814

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

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

Download in other formats:

Original Format