Context Navigation

watercloud.F @ 89

Last change on this file since 89 was 83, checked in by aslmd, 14 years ago

mars + LMD_MM_MARS : modifications pour cycle de l'eau, valeurs tunees JBM

used settings reached by JBM to obtain his PhD results

alb_surfice = 0.45 --- in physiq.F and meso_physiq.F
ccn_factor = 4.5 --- in watercloud.F
nuice_sed = 0.45 --- in callsedim.F

NB: this is supposed to be further refined in the future

File size: 11.5 KB

Rev	Line
[38]	1	SUBROUTINE watercloud(ngrid,nlay, ptimestep,
	2	& pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
	3	& pq,pdq,pdqcloud,pdqscloud,pdtcloud,
	4	& nq,naersize,tau,
	5	& ccn,rdust,rice,nuice)
	6	IMPLICIT NONE
	7
	8	c=======================================================================
	9	c Treatment of saturation of water vapor
	10	c
	11	c
	12	c Modif de zq si saturation dans l'atmosphere
	13	c si zq(ig,l)> zqsat(ig,l) -> zq(ig,l)=zqsat(ig,l)
	14	c Le test est effectue de bas en haut. L'eau condensee
	15	c (si saturation) est remise dans la couche en dessous.
	16	c L'eau condensee dans la couche du bas est deposee a la surface
	17	c
	18	c Modification: Franck Montmessin water ice scheme
	19	c Francois Forget : change nuclei density & outputs
	20	c Ehouarn Millour: sept.2008, tracers are now handled
	21	c by name (and not fixed index)
	22	c
	23	c=======================================================================
	24
	25	c-----------------------------------------------------------------------
	26	c declarations:
	27	c -------------
	28
	29	#include "dimensions.h"
	30	#include "dimphys.h"
	31	#include "comcstfi.h"
	32	#include "callkeys.h"
	33	#include "tracer.h"
	34	#include "comgeomfi.h"
	35
	36	c Inputs:
	37	c ------
	38
	39	INTEGER ngrid,nlay
	40	REAL ptimestep ! pas de temps physique (s)
	41	REAL pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
	42	REAL pplay(ngrid,nlay) ! pression au milieu des couches (Pa)
	43	REAL pdpsrf(ngrid) ! tendance surf pressure
	44	REAL pzlev(ngrid,nlay+1) ! altitude at layer boundaries
	45	REAL pzlay(ngrid,nlay) ! altitude at the middle of the layers
	46	REAL pt(ngrid,nlay) ! temperature at the middle of the layers (K)
	47	REAL pdt(ngrid,nlay) ! tendance temperature des autres param.
	48
	49	real pq(ngrid,nlay,nq) ! traceur (kg/kg)
	50	real pdq(ngrid,nlay,nq) ! tendance avant condensation (kg/kg.s-1)
	51
	52	integer nq ! nombre de traceurs
	53	integer naersize ! nombre de traceurs radiativement actifs (=naerkind)
	54	REAL tau(ngridmx,naersize)
	55	REAL ccn(ngridmx,nlayermx) ! Cloud condensation nuclei
	56	! (particules kg-1)
	57	real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)
	58
	59	c Outputs:
	60	c -------
	61
	62	real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation H2O(kg/kg.s-1)
	63	real pdqscloud(ngrid,nq) ! flux en surface (kg.m-2.s-1)
	64	REAL pdtcloud(ngrid,nlay) ! tendance temperature due
	65	! a la chaleur latente
	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
	71
	72	c local:
	73	c ------
	74
	75	REAL CBRT
	76	EXTERNAL CBRT
	77	INTEGER ig,iq,l
	78
	79
	80	REAL zq(ngridmx,nlayermx,nqmx) ! local value of tracers
	81	REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
	82	REAL zqsat(ngridmx,nlayermx) ! saturation
	83	REAL zt(ngridmx,nlayermx) ! local value of temperature
	84
	85	REAL masse (ngridmx,nlayermx)
	86	REAL epaisseur (ngridmx,nlayermx)
	87	REAL rfinal ! Ice crystal radius after condensation(m)
	88	REAL seq ! Equilibrium saturation ration (accounting for curvature effect)
	89	REAL dzq ! masse de glace echangee (kg/kg)
	90	REAL lw !Latent heat of sublimation (J.kg-1)
	91	REAL,PARAMETER :: To=273.15 ! reference temperature, T=273.15 K
	92
	93	REAL Ctot
	94	REAL*8 ph2o,satu
	95	REAL gr,Cste,up,dwn,newvap
	96
	97	LOGICAL,SAVE :: firstcall=.true.
	98	! To use more refined microphysics, set improved to .true.
	99	LOGICAL,PARAMETER :: improved=.true.
	100
	101	c Pour diagnostique :
	102	c ~~~~~~~~~~~~~~~~~
	103	c REAL icetot(ngridmx) ! Total mass of water ice (kg/m2)
	104	c REAL rave(ngridmx) ! Mean crystal radius in a column (m)
	105
	106	INTEGER i
	107
	108	! indexes of water vapour, water ice and dust tracers:
	109	INTEGER,SAVE :: i_h2o=0 ! water vapour
	110	INTEGER,SAVE :: i_ice=0 ! water ice
	111	CHARACTER(LEN=20) :: tracername ! to temporarly store text
	112
	113	c ** un petit test de coherence
	114	c --------------------------
	115
	116	IF (firstcall) THEN
	117	IF(ngrid.NE.ngridmx) THEN
	118	PRINT*,'STOP dans watercloud'
	119	PRINT*,'probleme de dimensions :'
	120	PRINT*,'ngrid =',ngrid
	121	PRINT*,'ngridmx =',ngridmx
	122	STOP
	123	ENDIF
	124
	125	if (nq.gt.nqmx) then
	126	write(,) 'stop in watercloud (nq.gt.nqmx)!'
	127	write(,) 'nq=',nq,' nqmx=',nqmx
	128	stop
	129	endif
	130
	131	i_h2o=igcm_h2o_vap
	132	i_ice=igcm_h2o_ice
	133
	134	write(,) "watercloud: i_h2o=",i_h2o
	135	write(,) " i_ice=",i_ice
	136
	137	firstcall=.false.
	138	ENDIF ! of IF (firstcall)
	139
	140
	141	c-----------------------------------------------------------------------
	142	c 1. initialisation
	143	c -----------------
	144
	145	c On "update" la valeur de q(nqmx) (water vapor) et temperature.
	146	c On effectue qqes calculs preliminaires sur les couches :
	147	c masse (kg.m-2), epaisseur(m).
	148
	149	do l=1,nlay
	150	do ig=1,ngrid
	151	zq(ig,l,i_h2o)=pq(ig,l,i_h2o)+pdq(ig,l,i_h2o)*ptimestep
	152	zq(ig,l,i_h2o)=max(zq(ig,l,i_h2o),1.E-30) ! FF 12/2004
	153	zq0(ig,l,i_h2o)=zq(ig,l,i_h2o)
	154	zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
	155	masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g
	156	epaisseur(ig,l)= pzlev(ig,l+1) - pzlev(ig,l)
	157
	158	zq(ig,l,i_ice)=pq(ig,l,i_ice)+pdq(ig,l,i_ice)*ptimestep
	159	zq(ig,l,i_ice)=max(zq(ig,l,i_ice),0.) ! FF 12/2004
	160	zq0(ig,l,i_ice)=zq(ig,l,i_ice)
	161
	162	c This typical profile is not used anymore; rdust is
	163	c set up in updatereffrad.F.
	164	c rdust(ig,l)= max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9)
	165	enddo
	166	enddo
	167
	168	do l=1,nlay
	169	do ig=1,ngrid
	170	c Calcul du rayon moyen des particules de glace.
	171	c Hypothese : Dans une couche, la glace presente se
	172	c repartit uniformement autour du nbre de poussieres
	173	c dont le rayon moyen est prescrit par rdust.
	174	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	175	rice(ig,l)=CBRT( ( zq(ig,l,i_ice)/rho_ice+
	176	& ccn(ig,l)(4./3.)pirdust(ig,l)*3. )
	177	& / (ccn(ig,l)4./3.pi) )
	178	rice(ig,l)=max(rice(ig,l),rdust(ig,l))
	179	c Effective variance of the size distribution
	180	nuice(ig,l)=nuice_ref
	181	enddo ! of do ig=1,ngrid
	182	enddo ! of dol=1,nlay
	183
	184	pdqscloud(1:ngrid,1:nq)=0
	185	pdqcloud(1:ngrid,1:nlay,1:nq)=0
	186	pdtcloud(1:ngrid,1:nlay)=0
	187
	188	c icetot(1:ngrid)=0
	189	c rave(1:ngrid)=0
	190
	191	c ----------------------------------------------
	192	c
	193	c
	194	c Rapport de melange a saturation dans la couche l : -------
	195	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	196
	197	call watersat(ngridmx*nlayermx,zt,pplay,zqsat)
	198
	199	c taux de condensation (kg/kg/s-1) dans les differentes couches
	200	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	201
	202	c Iceparty is not used anymore: water=>iceparty (JBM).
	203	c if(iceparty) then
	204
	205	do l=1,nlay
	206	do ig=1,ngrid
	207
	208	IF (improved) then
	209	c Improved microphysics scheme
	210	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	211
	212	Ctot = zq(ig,l,i_h2o) + zq(ig,l,i_ice)
	213	ph2o = zq(ig,l,i_h2o) * 44. / 18. * pplay(ig,l)
	214	satu = zq(ig,l,i_h2o) / zqsat(ig,l)
	215
	216	call growthrate(ptimestep,zt(ig,l),pplay(ig,l),
	217	& ph2o,ph2o/satu,seq,rice(ig,l),gr)
	218	Cste = ptimestep * 4. * pi * rice(ig,l)
	219	* * rho_ice * ccn(ig,l)
	220	up = zq(ig,l,i_h2o) + Cste * gr * seq
	221	dwn = 1. + Cste * gr / zqsat(ig,l)
	222	newvap = min(up/dwn,Ctot)
	223
	224	gr = gr * ( newvap/zqsat(ig,l) - seq )
	225	dzq = min( max( Cste * gr,-zq(ig,l,i_ice) )
	226	* , zq(ig,l,i_h2o) )
	227
	228	c Nucleation (sat ratio must be larger than a critical value)
	229	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	230	if (satu.gt.1.) then
	231	if (satu.le.1.4.and.zq(ig,l,i_ice).lt.1.e-8)
	232	* dzq = 0.
	233	endif
	234
	235	ELSE
	236	c Old version
	237	c ~~~~~~~~~~~
	238	if (zq(ig,l,i_h2o).ge.zqsat(ig,l))then ! Condensation
	239	dzq=zq(ig,l,i_h2o)-zqsat(ig,l)
	240	elseif(zq(ig,l,i_h2o).lt.zqsat(ig,l))then ! Sublimation
	241	dzq=-min(zqsat(ig,l)-zq(ig,l,i_h2o),zq(ig,l,i_ice))
	242	endif
	243
	244	ENDIF ! of IF (improved)
	245
	246	c Water Mass change
	247	c ~~~~~~~~~~~~~~~~~
	248	zq(ig,l,i_ice)=zq(ig,l,i_ice)+dzq
	249	zq(ig,l,i_h2o)=zq(ig,l,i_h2o)-dzq
	250
	251	rice(ig,l)=max( CBRT ( (zq(ig,l,i_ice)/rho_ice
	252	& +ccn(ig,l)(4./3.)pirdust(ig,l)*3.)
	253	& /(ccn(ig,l)4./3.pi) ), rdust(ig,l))
	254
	255	enddo ! of do ig=1,ngrid
	256	enddo ! of do l=1,nlay
	257
	258	c The following part have been commented because iceparty
	259	c is not used anymore: water=>iceparty (JBM).
	260
	261	c else ! if not iceparty
	262
	263	c Saturation couche nlay a 2 :
	264	c ~~~~~~~~~~~~~~~~~~~~~~~~~~
	265	c do l=nlay,2, -1
	266	c do ig=1,ngrid
	267	c if (zq(ig,l,i_h2o).gt.zqsat(ig,l))then
	268	c zq(ig,l-1,i_h2o)= zq(ig,l-1,i_h2o)+
	269	c & (zq(ig,l,i_h2o)-zqsat(ig,l))
	270	c & *(pplev(ig,l)-pplev(ig,l+1))/(pplev(ig,l-1)-pplev(ig,l))
	271	c zq(ig,l,i_h2o)=zqsat(ig,l)
	272	c endif
	273	c enddo
	274	c enddo
	275
	276	c Saturation couche l=1 si pas iceparty
	277	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	278	c do ig=1,ngridmx
	279	c if (zq(ig,1,i_h2o).gt.zqsat(ig,1))then
	280	c pdqscloud(ig,i_ice)=(zq(ig,1,i_h2o)-zqsat(ig,1))
	281	c & (pplev(ig,1)-pplev(ig,2))/(gptimestep)
	282	c zq(ig,1,i_h2o)=zqsat(ig,1)
	283	c endif
	284	c enddo
	285
	286	c endif ! of if (iceparty)
	287
	288	c Tendance finale
	289	c ~~~~~~~~~~~~~~~
	290	do l=1, nlay
	291	do ig=1,ngridmx
	292	pdqcloud(ig,l,i_h2o)=(zq(ig,l,i_h2o)
	293	& -zq0(ig,l,i_h2o))/ptimestep
	294	pdqcloud(ig,l,i_ice) =
	295	& (zq(ig,l,i_ice) - zq0(ig,l,i_ice))/ptimestep
	296	lw=(2834.3-0.28(zt(ig,l)-To)-0.004(zt(ig,l)-To)*2)1.e+3
	297	pdtcloud(ig,l)=-pdqcloud(ig,l,i_h2o)*lw/cpp
	298	end do
	299	end do
	300
	301	c A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
	302	c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	303	c Then that should not affect the ice particle radius
	304	do ig=1,ngridmx
	305	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,2)))then
	306	if(pdpsrf(ig)ptimestep.gt.0.9(pplev(ig,1)-pplev(ig,3)))
	307	& rice(ig,2)=rice(ig,3)
	308	rice(ig,1)=rice(ig,2)
	309	end if
	310	end do
	311
	312	c**************************************************
	313	c Output
	314	c**************************************************
	315	! NB: for diagnostics use zq(), the updated value of tracers
	316
	317	c do ig=1,ngridmx
	318	c do l=1 ,nlay
	319	c masse de glace d'eau dans la couche l
	320	c icetot(ig)=icetot(ig)+masse(ig,l)*zq(ig,l,i_ice)
	321	c rayon moyen des cristaux dans la colonne ig
	322	c rave(ig)=rave(ig)+masse(ig,l)zq(ig,l,i_ice)rice(ig,l)
	323	c enddo
	324	c rave(ig)=rave(ig)/max(icetot(ig),1.e-30)
	325	c if (icetot(ig)*1000.lt.0.01) rave(ig)=0.
	326	c enddo ! (ngridmx)
	327	c**************************************************
	328
	329	RETURN
	330	END
	331

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

Download in other formats:

Original Format