Context Navigation

watercloud_mod.F @ 2932

Last change on this file since 2932 was 2932, checked in by romain.vande, 21 months ago

Mars PCM:
Add a new routine to write output called write_output.
It needs to be imported (for example : use write_output_mod, only: write_output)
Then, it requires only 4 arguments : the name of the variable, its title, its units and the variable itself.
It detects the dimension of the variable and decide to output either in diagfi or diagsoil.
It is also compatible with XIOS (xml file needs to be adapted)
Writediagfi and writediagsoil routines are still available in case.
RV

File size: 28.5 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
[2932]	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=======================================================================
[2932]	716	call write_output("pdqice2","pdqcloudice apres microphysique"
	717	& ,"kg/kg.s-1",pdqcloud(:,:,igcm_h2o_ice))
	718	call write_output("pdqvap2","pdqcloudvap apres microphysique"
	719	& ,"kg/kg.s-1",pdqcloud(:,:,
[1711]	720	& igcm_h2o_vap))
[2312]	721	if (hdo) then
[2932]	722	call write_output("pdqiceD","pdqiceD apres microphysique"
	723	& ,"kg/kg.s-1",pdqcloud(:,:,igcm_hdo_ice))
	724	call write_output("pdqvapD","pdqvapD apres microphysique"
	725	& ,"kg/kg.s-1",pdqcloud(:,:,
[2312]	726	& igcm_hdo_vap))
	727	endif
[2932]	728	call write_output("pdqccn2","pdqcloudccn apres microphysique"
	729	& ,"kg/kg.s-1",pdqcloud(:,:,
[1711]	730	& igcm_ccn_mass))
[2932]	731	call write_output("pdqccnN2","pdqcloudccnN apres "//
	732	& "microphysique","nb/kg.s-1",pdqcloud(:,:,
[1711]	733	& igcm_ccn_number))
[2932]	734	call write_output("pdqdust2", "pdqclouddust apres "//
	735	& "microphysique","kg/kg.s-1",pdqcloud(:,:,
[1711]	736	& igcm_dust_mass))
[2932]	737	call write_output("pdqdustN2", "pdqclouddustN apres "//
	738	& "microphysique","nb/kg.s-1",pdqcloud(:,:,
[1711]	739	& igcm_dust_number))
[633]	740	c=======================================================================
[1758]	741	#endif
[633]	742
[1711]	743	END SUBROUTINE watercloud
	744
[2162]	745	subroutine ini_watercloud_mod(ngrid,nlayer,nq)
	746	implicit none
	747
	748	integer,intent(in) :: ngrid ! number of atmospheric columns
	749	integer,intent(in) :: nlayer ! number of atmospheric layers
	750	integer,intent(in) :: nq ! number of tracers
	751
	752	allocate(zdqcloud(ngrid,nlayer,nq))
	753	zdqcloud(:,:,:)=0
	754	allocate(zdqscloud(ngrid,nq))
	755	zdqscloud(:,:)=0
	756
	757	end subroutine ini_watercloud_mod
	758
	759
	760	subroutine end_watercloud_mod
	761	implicit none
	762
	763	if (allocated(zdqcloud)) deallocate(zdqcloud)
	764	if (allocated(zdqscloud)) deallocate(zdqscloud)
	765
	766	end subroutine end_watercloud_mod
	767
[1711]	768	END MODULE watercloud_mod

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

Download in other formats:

Original Format