Context Navigation

watercloud.F @ 522

Last change on this file since 522 was 522, checked in by tnavarro, 13 years ago
forgot the most important...
File size: 16.5 KB

Rev	Line
[522]	1	SUBROUTINE watercloud(ngrid,nlay,ptimestep,
[38]	2	& pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
	3	& pq,pdq,pdqcloud,pdqscloud,pdtcloud,
[358]	4	& nq,tau,tauscaling,rdust,rice,nuice,
	5	& rsedcloud,rhocloud)
[522]	6	! to use 'getin'
	7	USE ioipsl_getincom
[38]	8	IMPLICIT NONE
	9
	10	c=======================================================================
[358]	11	c Water-ice cloud formation
	12	c
	13	c Includes two different schemes:
	14	c - A simplified scheme (see simpleclouds.F)
	15	c - An improved microphysical scheme (see improvedclouds.F)
[38]	16	c
[522]	17	c There is a time loop specific to cloud formation and sedimentation
	18	c due to timescales smaller than the GCM integration timestep.
	19	c
[358]	20	c Authors: Franck Montmessin, Francois Forget, Ehouarn Millour,
[522]	21	c J.-B. Madeleine, Thomas Navarro
[38]	22	c
[522]	23	c 2004 - 2012
[38]	24	c=======================================================================
	25
	26	c-----------------------------------------------------------------------
	27	c declarations:
	28	c -------------
	29
	30	#include "dimensions.h"
	31	#include "dimphys.h"
	32	#include "comcstfi.h"
	33	#include "callkeys.h"
	34	#include "tracer.h"
	35	#include "comgeomfi.h"
[358]	36	#include "dimradmars.h"
[38]	37
	38	c Inputs:
	39	c ------
	40
	41	INTEGER ngrid,nlay
[522]	42	INTEGER nq ! nombre de traceurs
[38]	43	REAL ptimestep ! pas de temps physique (s)
	44	REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
	45	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
[522]	46	REAL pdpsrf(ngrid) ! tendence surf pressure
[38]	47	REAL pzlev(ngrid,nlay+1) ! altitude at layer boundaries
	48	REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers
	49	REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K)
[522]	50	REAL pdt(ngrid,nlay) ! tendence temperature des autres param.
[38]	51
	52	real pq(ngrid,nlay,nq) ! traceur (kg/kg)
[522]	53	real pdq(ngrid,nlay,nq) ! tendence avant condensation (kg/kg.s-1)
[38]	54
[522]	55	REAL tau(ngridmx,naerkind) ! Column dust optical depth at each point
	56	REAL tauscaling(ngridmx) ! Convertion factor for dust amount
	57	real rdust(ngridmx,nlay) ! Dust geometric mean radius (m)
[38]	58
	59	c Outputs:
	60	c -------
	61
[522]	62	real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
[38]	63	real pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1)
[522]	64	REAL pdtcloud(ngrid,nlay) ! tendence temperature due
	65	! a la chaleur latente
[38]	66
	67	REAL rice(ngrid,nlay) ! Ice mass mean radius (m)
	68	! (r_c in montmessin_2004)
	69	REAL nuice(ngrid,nlay) ! Estimated effective variance
	70	! of the size distribution
[522]	71	real rsedcloud(ngridmx,nlay) ! Cloud sedimentation radius
	72	real rhocloud(ngridmx,nlay) ! Cloud density (kg.m-3)
[38]	73
	74	c local:
	75	c ------
	76
[522]	77	! for sedimentation
	78	REAL zq(ngridmx,nlay,nqmx) ! local value of tracers
	79	real masse (ngridmx,nlay) ! Layer mass (kg.m-2)
	80	real epaisseur (ngridmx,nlay) ! Layer thickness (m)
	81	real wq(ngridmx,nlay+1) ! displaced tracer mass (kg.m-2)
	82
	83	! for ice radius computation
	84	REAL Mo,No
	85	REAl ccntyp
	86	real beta ! correction for the shape of the ice particles (cf. newsedim)
	87	save beta
	88
	89	! for time loop
	90	INTEGER microstep ! time subsampling step variable
	91	INTEGER imicro ! time subsampling for coupled water microphysics & sedimentation
	92	SAVE imicro
	93	REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
	94	SAVE microtimestep
	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
	99
	100	! global tendency (clouds+sedim+physics)
	101	REAL subpdq(ngrid,nlay,nq) ! cf. pdqcloud
	102	REAL subpdt(ngrid,nlay) ! cf. pdtcloud
	103
	104	REAL CBRT
	105	EXTERNAL CBRT
	106
	107
	108	INTEGER iq,ig,l
[38]	109	LOGICAL,SAVE :: firstcall=.true.
	110
	111	c ** un petit test de coherence
	112	c --------------------------
	113
	114	IF (firstcall) THEN
	115	IF(ngrid.NE.ngridmx) THEN
	116	PRINT*,'STOP dans watercloud'
	117	PRINT*,'probleme de dimensions :'
	118	PRINT*,'ngrid =',ngrid
	119	PRINT*,'ngridmx =',ngridmx
	120	STOP
	121	ENDIF
	122
	123	if (nq.gt.nqmx) then
	124	write(,) 'stop in watercloud (nq.gt.nqmx)!'
	125	write(,) 'nq=',nq,' nqmx=',nqmx
	126	stop
	127	endif
	128
[358]	129	write(,) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
	130	write(,) " igcm_h2o_ice=",igcm_h2o_ice
[522]	131
[38]	132
[522]	133	write(,) "correction for the shape of the ice particles ?"
	134	beta=0.75 ! default value
	135	call getin("ice_shape",beta)
	136	write(,) " ice_shape = ",beta
	137
	138	write(,) "time subsampling for microphysic ?"
	139	imicro = 1
	140	call getin("imicro",imicro)
	141	write(,)"imicro = ",imicro
	142
	143	microtimestep = ptimestep/float(imicro)
	144	write(,)"Physical timestep is",ptimestep
	145	write(,)"Microphysics timestep is",microtimestep
	146
[38]	147	firstcall=.false.
	148	ENDIF ! of IF (firstcall)
[522]	149
	150	c-----Initialization
	151	subpdq(1:ngrid,1:nlay,1:nq) = 0
	152	subpdt(1:ngrid,1:nlay) = 0
	153	pdqscloud(1:ngrid,1:nq) = 0
	154	zq(1:ngrid,1:nlay,1:nq) = 0
[38]	155
[522]	156	c-----Computing the different layer properties for clouds sedimentation
	157	do l=1,nlay
	158	do ig=1, ngrid
	159	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
	160	epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
	161	enddo
	162	enddo
	163
	164	c------------------------------------------------------------------
	165	c Time subsampling for coupled microphysics and sedimentation
	166	c------------------------------------------------------------------
	167	DO microstep=1,imicro
	168
	169	c-------------------------------------------------------------------
	170	c 1. Tendencies:
	171	c------------------
	172
	173
	174	c------ Temperature tendency subpdt
	175	! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
	176	! If imicro=1 subpdt is the same as pdt
	177	DO l=1,nlay
	178	DO ig=1,ngrid
	179	subpdt(ig,l) = subpdt(ig,l)
	180	& + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
	181	ENDDO
	182	ENDDO
	183	c------ Traceurs tendencies subpdq
	184	c------ At each micro timestep we add pdq in order to have a stepped entry
	185	IF (microphys) THEN
	186	DO l=1,nlay
	187	DO ig=1,ngrid
	188	subpdq(ig,l,igcm_dust_mass) =
	189	& subpdq(ig,l,igcm_dust_mass)
	190	& + pdq(ig,l,igcm_dust_mass)
	191	subpdq(ig,l,igcm_dust_number) =
	192	& subpdq(ig,l,igcm_dust_number)
	193	& + pdq(ig,l,igcm_dust_number)
	194	subpdq(ig,l,igcm_ccn_mass) =
	195	& subpdq(ig,l,igcm_ccn_mass)
	196	& + pdq(ig,l,igcm_ccn_mass)
	197	subpdq(ig,l,igcm_ccn_number) =
	198	& subpdq(ig,l,igcm_ccn_number)
	199	& + pdq(ig,l,igcm_ccn_number)
	200	ENDDO
	201	ENDDO
	202	ENDIF
	203	DO l=1,nlay
	204	DO ig=1,ngrid
	205	subpdq(ig,l,igcm_h2o_ice) =
	206	& subpdq(ig,l,igcm_h2o_ice)
	207	& + pdq(ig,l,igcm_h2o_ice)
	208	subpdq(ig,l,igcm_h2o_vap) =
	209	& subpdq(ig,l,igcm_h2o_vap)
	210	& + pdq(ig,l,igcm_h2o_vap)
	211	ENDDO
	212	ENDDO
	213
	214
	215	c-------------------------------------------------------------------
	216	c 2. Main call to the different cloud schemes:
	217	c------------------------------------------------
	218	IF (microphys) THEN
	219	CALL improvedclouds(ngrid,nlay,microtimestep,
	220	& pplev,pplay,pzlev,pt,subpdt,
	221	& pq,subpdq,subpdqcloud,subpdtcloud,
[358]	222	& nq,tauscaling,rdust,rice,nuice,
	223	& rsedcloud,rhocloud)
[522]	224	ELSE
	225	CALL simpleclouds(ngrid,nlay,microtimestep,
	226	& pplev,pplay,pzlev,pzlay,pt,subpdt,
	227	& pq,subpdq,subpdqcloud,subpdtcloud,
[358]	228	& nq,tau,rice,nuice,rsedcloud)
[522]	229	ENDIF
	230
	231	c--------------------------------------------------------------------
	232	c 3. CCN & ice sedimentation:
	233	c--------------------------------
	234	! Sedimentation is done here for water ice and its CCN only
	235	! callsedim in physics is done for dust (and others species if any)
	236
	237	DO l=1,nlay
	238	DO ig=1,ngrid
	239	subpdt(ig,l) =
	240	& subpdt(ig,l) + subpdtcloud(ig,l)
	241	ENDDO
	242	ENDDO
	243
	244	c------- water ice update before sedimentation (radius is done in the cloud scheme routine)
	245	DO l=1,nlay
	246	DO ig=1, ngrid
	247	zq(ig,l,igcm_h2o_ice)= max(1e-30,
	248	& pq(ig,l,igcm_h2o_ice) + (subpdq(ig,l,igcm_h2o_ice)
	249	& + subpdqcloud(ig,l,igcm_h2o_ice))*microtimestep)
	250	! zq(ig,l,igcm_h2o_vap)= max(1e-30,
	251	! & pq(ig,l,igcm_h2o_vap) + (subpdq(ig,l,igcm_h2o_vap)
	252	! & + subpdqcloud(ig,l,igcm_h2o_vap))*microtimestep)
	253	ENDDO
	254	ENDDO
	255
	256
	257	c------- CCN update before sedimentation
	258	IF (microphys) THEN
	259	DO l=1,nlay
	260	DO ig=1,ngrid
	261	zq(ig,l,igcm_ccn_number)=
	262	& pq(ig,l,igcm_ccn_number) + (subpdq(ig,l,igcm_ccn_number)
	263	& + subpdqcloud(ig,l,igcm_ccn_number))*microtimestep
	264	zq(ig,l,igcm_ccn_number)= max(1e-30,
	265	& zq(ig,l,igcm_ccn_number))!*tauscaling(ig)) ! OU pas tauscaling ?
	266	zq(ig,l,igcm_ccn_mass)=
	267	& pq(ig,l,igcm_ccn_mass) + (subpdq(ig,l,igcm_ccn_mass)
	268	& + subpdqcloud(ig,l,igcm_ccn_mass))*microtimestep
	269	zq(ig,l,igcm_ccn_mass)=max(1e-30,
	270	& zq(ig,l,igcm_ccn_mass))!*tauscaling(ig)) ! OU pas tauscaling ?
	271	ENDDO
	272	ENDDO
	273	ENDIF
	274
[38]	275
[522]	276
	277	IF (microphys) THEN
	278
	279	c------- CCN number sedimentation
	280	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
	281	& microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
	282	& rhocloud,zq(1,1,igcm_ccn_number),wq,beta)
	283	do ig=1,ngrid
	284	! matters if one would like to know ccn surface deposition
	285	pdqscloud(ig,igcm_ccn_number)=
	286	& pdqscloud(ig,igcm_ccn_number) + wq(ig,1)
	287	enddo
	288
	289	c------- CCN mass sedimentation
	290	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
	291	& microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
	292	& rhocloud,zq(1,1,igcm_ccn_mass),wq,beta)
	293	do ig=1,ngrid
	294	! matters if one would like to know ccn surface deposition
	295	pdqscloud(ig,igcm_ccn_mass)=
	296	& pdqscloud(ig,igcm_ccn_mass) + wq(ig,1)
	297	enddo
	298
	299	c------- Water ice sedimentation -- improved microphys
	300	call newsedim(ngrid,nlay,ngridnlay,ngridnlay,
	301	& microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
	302	& rhocloud,zq(1,1,igcm_h2o_ice),wq,beta)
	303	ELSE
	304
	305	c------- Water ice sedimentation -- simple microphys
	306	call newsedim(ngrid,nlay,ngrid*nlay,1,
	307	& microtimestep,pplev,masse,epaisseur,pt,rsedcloud,
	308	& rho_q(igcm_h2o_ice),zq(1,1,igcm_h2o_ice),wq,beta)
	309
	310	ENDIF
	311
	312
	313	c------- Surface ice tendency update
	314	DO ig=1,ngrid
	315	pdqscloud(ig,igcm_h2o_ice)=
	316	& pdqscloud(ig,igcm_h2o_ice) + wq(ig,1)
	317	ENDDO
	318
	319
	320	c-------------------------------------------------------------------
	321	c 5. Updating tendencies after sedimentation:
	322	c-----------------------------------------------
	323
	324	DO l=1,nlay
	325	DO ig=1,ngrid
	326
	327	subpdq(ig,l,igcm_h2o_ice) =
	328	& (zq(ig,l,igcm_h2o_ice) - pq(ig,l,igcm_h2o_ice))
	329	& /microtimestep
	330
	331
	332	subpdq(ig,l,igcm_h2o_vap)=subpdq(ig,l,igcm_h2o_vap)
	333	& +subpdqcloud(ig,l,igcm_h2o_vap)
	334
	335	ENDDO
	336	ENDDO
	337
	338
	339	IF (microphys) then
	340	DO l=1,nlay
	341	DO ig=1,ngrid
	342	subpdq(ig,l,igcm_ccn_number)=(zq(ig,l,igcm_ccn_number)
	343	& -pq(ig,l,igcm_ccn_number))/microtimestep
	344	subpdq(ig,l,igcm_ccn_mass)=(zq(ig,l,igcm_ccn_mass)
	345	& -pq(ig,l,igcm_ccn_mass))/microtimestep
	346	ENDDO
	347	ENDDO
	348	ENDIF
	349
	350
	351
	352	ENDDO ! of DO microstep=1,imicro
	353
	354	c-------------------------------------------------------------------
	355	c 6. Compute final tendencies after time loop:
	356	c------------------------------------------------
	357
	358	c------ Whole temperature tendency pdtcloud
	359	DO l=1,nlay
	360	DO ig=1,ngrid
	361	pdtcloud(ig,l) =
	362	& subpdt(ig,l)/imicro-pdt(ig,l)
	363	ENDDO
	364	ENDDO
	365	c------ Traceurs tendencies pdqcloud
	366	DO iq=1,nq
	367	DO l=1,nlay
	368	DO ig=1,ngrid
	369	pdqcloud(ig,l,iq) = subpdq(ig,l,iq)/imicro
	370	& - pdq(ig,l,iq)
	371	ENDDO
	372	ENDDO
	373	ENDDO
	374	c------ Traceurs surface tendencies pdqscloud
	375	DO iq=1,nq
	376	DO ig=1,ngrid
	377	pdqscloud(ig,iq) =
	378	& pdqscloud(ig,iq)/ptimestep
	379	ENDDO
	380	ENDDO
	381
	382
	383
	384	c------Update the ice particle size "rice" for output or photochemistry.
	385	c------It is not used again in the water cycle.
	386	IF(scavenging) THEN
	387	DO l=1, nlay
	388	DO ig=1,ngrid
	389	Mo = pq(ig,l,igcm_h2o_ice)
	390	& + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep
	391	& + (pq(ig,l,igcm_ccn_mass)
	392	& + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
	393	& *tauscaling(ig) + 1.e-30
	394	No = (pq(ig,l,igcm_ccn_number)
	395	& + pdqcloud(ig,l,igcm_ccn_number)*ptimestep)
	396	& *tauscaling(ig) + 1.e-30
	397	rhocloud(ig,l) = (pq(ig,l,igcm_h2o_ice) +
	398	& pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)
	399	& / Mo * rho_ice
	400	& + (pq(ig,l,igcm_ccn_mass)
	401	& + pdqcloud(ig,l,igcm_ccn_mass)*ptimestep)
	402	& tauscaling(ig)/ Mo rho_dust
	403	rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
	404	rice(ig,l)=(Mo / No * 0.75 / pi / rhocloud(ig,l))**(1./3.)
	405	if ((Mo.lt.1.e-15) .or. (No.le.50)) rice(ig,l) = 1.e-8
	406	ENDDO
	407	ENDDO
	408	ELSE
	409	DO l=1,nlay
	410	DO ig=1,ngrid
	411	ccntyp =
	412	& 1.3e+8max(tau(ig,1),0.001)/0.1exp(-pzlay(ig,l)/10000.)
	413	ccntyp = ccntyp /ccn_factor
	414	rice(ig,l)=max( CBRT ( ((pq(ig,l,igcm_h2o_ice)
	415	& + pdqcloud(ig,l,igcm_h2o_ice)*ptimestep)/rho_ice
	416	& +ccntyp(4./3.)pirdust(ig,l)*3.)
	417	& /(ccntyp4./3.pi) ), rdust(ig,l))
	418	ENDDO
	419	ENDDO
	420	ENDIF ! of IF(scavenging)
	421
	422	!--------------------------------------------------------------
	423	!--------------------------------------------------------------
	424
	425
[358]	426	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
	427	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	428	c Then that should not affect the ice particle radius
	429	do ig=1,ngridmx
	430	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
	431	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
	432	& rice(ig,2)=rice(ig,3)
	433	rice(ig,1)=rice(ig,2)
	434	end if
	435	end do
[38]	436
[522]	437	c=======================================================================
	438
	439	!!!!!!!!!! FOR PHOTOCHEMISTRY, REIMPLEMENT output of surfdust/surfice
	440	!! if (photochem) then
	441	!!c computation of dust and ice surface area (micron2/cm3)
	442	!!c for heterogeneous chemistry
	443	!!
	444	!! do l = 1,nlay
	445	!! do ig = 1,ngrid
	446	!!c
	447	!!c npart: number density of ccn in #/cm3
	448	!!c
	449	!! npart(ig,l) = 1.e-6*ccn(ig,l)
	450	!! $ *masse(ig,l)/epaisseur(ig,l)
	451	!!c
	452	!!c dust and ice surface area
	453	!!c
	454	!! surfdust(ig,l) = npart(ig,l)4.pi1.e12rdust(ig,l)**2
	455	!!c
	456	!! if (rice(ig,l) .ge. rdust(ig,l)) then
	457	!! surfice(ig,l) = npart(ig,l)4.pi1.e12rice(ig,l)**2
	458	!! surfdust(ig,l) = 0.
	459	!! else
	460	!! surfice(ig,l) = 0.
	461	!! end if
	462	!! end do ! of do ig=1,ngrid
	463	!! end do ! of do l=1,nlay
	464	!! end if ! of photochem
	465
[38]	466	RETURN
	467	END
	468

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

Download in other formats:

Original Format