Context Navigation

rain.F90 @ 1350

Last change on this file since 1350 was 1315, checked in by milmd, 12 years ago
LMDZ.GENERIC. OpenMP directives added in generic physic. When running in pure OpenMP or hybrid OpenMP/MPI, may have some bugs with condense_cloud and wstats routines.
File size: 14.1 KB

Rev	Line
[1308]	1	subroutine rain(ngrid,nlayer,nq,ptimestep,pplev,pplay,t,pdt,pq,pdq,d_t,dqrain,dqsrain,dqssnow,rneb)
[135]	2
[253]	3
[1315]	4	! to use 'getin'
	5	! use ioipsl_getincom
	6	use ioipsl_getincom_p
[728]	7	use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds, RLVTT, RCPD, RCPV, RV, RVTMP2,Psat_water,Tsat_water,rhowater
	8	use radii_mod, only: h2o_cloudrad
[1283]	9	USE tracer_h, only: igcm_h2o_vap, igcm_h2o_ice
[135]	10	implicit none
	11
	12	!==================================================================
	13	!
	14	! Purpose
	15	! -------
	16	! Calculates H2O precipitation using simplified microphysics.
	17	!
	18	! Authors
	19	! -------
	20	! Adapted from the LMDTERRE code by R. Wordsworth (2009)
[728]	21	! Added rain vaporization in case of T>Tsat
[135]	22	! Original author Z. X. Li (1993)
	23	!
	24	!==================================================================
	25
[1308]	26	! include "dimensions.h"
	27	! include "dimphys.h"
[1283]	28	include "comcstfi.h"
	29	include "callkeys.h"
[135]	30
[858]	31	! Arguments
[1308]	32	integer,intent(in) :: ngrid ! number of atmospheric columns
	33	integer,intent(in) :: nlayer ! number of atmospheric layers
[858]	34	integer,intent(in) :: nq ! number of tracers
	35	real,intent(in) :: ptimestep ! time interval
[1308]	36	real,intent(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
	37	real,intent(in) :: pplay(ngrid,nlayer) ! mid-layer pressure (Pa)
	38	real,intent(in) :: t(ngrid,nlayer) ! input temperature (K)
	39	real,intent(in) :: pdt(ngrid,nlayer) ! input tendency on temperature (K/s)
	40	real,intent(in) :: pq(ngrid,nlayer,nq) ! tracers (kg/kg)
	41	real,intent(in) :: pdq(ngrid,nlayer,nq) ! input tendency on tracers
	42	real,intent(out) :: d_t(ngrid,nlayer) ! temperature tendency (K/s)
	43	real,intent(out) :: dqrain(ngrid,nlayer,nq) ! tendency of H2O precipitation (kg/kg.s-1)
[858]	44	real,intent(out) :: dqsrain(ngrid) ! rain flux at the surface (kg.m-2.s-1)
	45	real,intent(out) :: dqssnow(ngrid) ! snow flux at the surface (kg.m-2.s-1)
[1308]	46	real,intent(in) :: rneb(ngrid,nlayer) ! cloud fraction
[787]	47
[1308]	48	REAL zt(ngrid,nlayer) ! working temperature (K)
	49	REAL ql(ngrid,nlayer) ! liquid water (Kg/Kg)
	50	REAL q(ngrid,nlayer) ! specific humidity (Kg/Kg)
	51	REAL d_q(ngrid,nlayer) ! water vapor increment
	52	REAL d_ql(ngrid,nlayer) ! liquid water / ice increment
[135]	53
	54	! Subroutine options
[858]	55	REAL,PARAMETER :: seuil_neb=0.001 ! Nebulosity threshold
[135]	56
[728]	57	INTEGER,save :: precip_scheme ! id number for precipitaion scheme
	58	! for simple scheme (precip_scheme=1)
	59	REAL,SAVE :: rainthreshold ! Precipitation threshold in simple scheme
	60	! for sundquist scheme (precip_scheme=2-3)
	61	REAL,SAVE :: cloud_sat ! Precipitation threshold in non simple scheme
	62	REAL,SAVE :: precip_timescale ! Precipitation timescale
	63	! for Boucher scheme (precip_scheme=4)
	64	REAL,SAVE :: Cboucher ! Precipitation constant in Boucher 95 scheme
	65	REAL,PARAMETER :: Kboucher=1.19E8
	66	REAL,SAVE :: c1
[1315]	67	!$OMP THREADPRIVATE(precip_scheme,rainthreshold,cloud_sat,precip_timescale,Cboucher,c1)
[135]	68
[858]	69	INTEGER,PARAMETER :: ninter=5
[135]	70
[1016]	71	logical,save :: evap_prec ! Does the rain evaporate?
[1315]	72	!$OMP THREADPRIVATE(evap_prec)
[135]	73
	74	! for simple scheme
[858]	75	real,parameter :: t_crit=218.0
[135]	76	real lconvert
	77
	78	! Local variables
	79	INTEGER i, k, n
[1308]	80	REAL zqs(ngrid,nlayer),Tsat(ngrid,nlayer), zdelta, zcor
[787]	81	REAL zrfl(ngrid), zrfln(ngrid), zqev, zqevt
[135]	82
[787]	83	REAL zoliq(ngrid)
	84	REAL zdz(ngrid),zrho(ngrid),ztot(ngrid), zrhol(ngrid)
	85	REAL zchau(ngrid),zfroi(ngrid),zfrac(ngrid),zneb(ngrid)
[135]	86
[1308]	87	real reffh2oliq(ngrid,nlayer),reffh2oice(ngrid,nlayer)
[728]	88
	89	real ttemp, ptemp, psat_tmp
[1308]	90	real tnext(ngrid,nlayer)
[135]	91
[1308]	92	real l2c(ngrid,nlayer)
[253]	93	real dWtot
[135]	94
[253]	95
[135]	96	! Indices of water vapour and water ice tracers
	97	INTEGER, SAVE :: i_vap=0 ! water vapour
	98	INTEGER, SAVE :: i_ice=0 ! water ice
[1315]	99	!$OMP THREADPRIVATE(i_vap,i_ice)
[135]	100
[858]	101	LOGICAL,SAVE :: firstcall=.true.
[1315]	102	!$OMP THREADPRIVATE(firstcall)
[135]	103
	104	! Online functions
[731]	105	REAL fallv, fall2v, zzz ! falling speed of ice crystals
[135]	106	fallv (zzz) = 3.29 * ((zzz)**0.16)
[731]	107	fall2v (zzz) =10.6 * ((zzz)**0.31) !for use with radii
[135]	108
	109
	110	IF (firstcall) THEN
	111
	112	i_vap=igcm_h2o_vap
	113	i_ice=igcm_h2o_ice
	114
	115	write(,) "rain: i_ice=",i_ice
	116	write(,) " i_vap=",i_vap
	117
	118	PRINT*, 'in rain.F, ninter=', ninter
	119	PRINT*, 'in rain.F, evap_prec=', evap_prec
	120
[728]	121	write(,) "Precipitation scheme to use?"
	122	precip_scheme=1 ! default value
[1315]	123	call getin_p("precip_scheme",precip_scheme)
[728]	124	write(,) " precip_scheme = ",precip_scheme
	125
	126	if (precip_scheme.eq.1) then
	127	write(,) "rainthreshold in simple scheme?"
	128	rainthreshold=0. ! default value
[1315]	129	call getin_p("rainthreshold",rainthreshold)
[728]	130	write(,) " rainthreshold = ",rainthreshold
	131
	132	else if (precip_scheme.eq.2.or.precip_scheme.eq.3) then
	133	write(,) "cloud water saturation level in non simple scheme?"
	134	cloud_sat=2.6e-4 ! default value
[1315]	135	call getin_p("cloud_sat",cloud_sat)
[728]	136	write(,) " cloud_sat = ",cloud_sat
	137	write(,) "precipitation timescale in non simple scheme?"
	138	precip_timescale=3600. ! default value
[1315]	139	call getin_p("precip_timescale",precip_timescale)
[728]	140	write(,) " precip_timescale = ",precip_timescale
	141
	142	else if (precip_scheme.eq.4) then
	143	write(,) "multiplicative constant in Boucher 95 precip scheme"
	144	Cboucher=1. ! default value
[1315]	145	call getin_p("Cboucher",Cboucher)
[728]	146	write(,) " Cboucher = ",Cboucher
	147	c1=1.001.097/rhowaterCboucher*Kboucher
	148
	149	endif
	150
[1016]	151	write(,) "re-evaporate precipitations?"
	152	evap_prec=.true. ! default value
[1315]	153	call getin_p("evap_prec",evap_prec)
[1016]	154	write(,) " evap_prec = ",evap_prec
	155
[135]	156	firstcall = .false.
[1283]	157	ENDIF ! of IF (firstcall)
[135]	158
	159	! GCM -----> subroutine variables
[1308]	160	DO k = 1, nlayer
[787]	161	DO i = 1, ngrid
[135]	162
[253]	163	zt(i,k) = t(i,k)+pdt(i,k)*ptimestep ! a big fat bug was here
	164	q(i,k) = pq(i,k,i_vap)+pdq(i,k,i_vap)*ptimestep
	165	ql(i,k) = pq(i,k,i_ice)+pdq(i,k,i_ice)*ptimestep
[135]	166
[253]	167	!q(i,k) = pq(i,k,i_vap)!+pdq(i,k,i_vap)
	168	!ql(i,k) = pq(i,k,i_ice)!+pdq(i,k,i_ice)
	169
[135]	170	if(q(i,k).lt.0.)then ! if this is not done, we don't conserve water
	171	q(i,k)=0.
	172	endif
	173	if(ql(i,k).lt.0.)then
	174	ql(i,k)=0.
	175	endif
	176
	177	ENDDO
	178	ENDDO
	179
	180	! Initialise the outputs
[1308]	181	d_t(1:ngrid,1:nlayer) = 0.0
	182	d_q(1:ngrid,1:nlayer) = 0.0
	183	d_ql(1:ngrid,1:nlayer) = 0.0
[1283]	184	zrfl(1:ngrid) = 0.0
	185	zrfln(1:ngrid) = 0.0
[135]	186
	187	! calculate saturation mixing ratio
[1308]	188	DO k = 1, nlayer
[787]	189	DO i = 1, ngrid
[253]	190	ttemp = zt(i,k)
[135]	191	ptemp = pplay(i,k)
[728]	192	! call watersat(ttemp,ptemp,zqs(i,k))
	193	call Psat_water(ttemp,ptemp,psat_tmp,zqs(i,k))
	194	call Tsat_water(ptemp,Tsat(i,k))
[135]	195	ENDDO
	196	ENDDO
	197
[253]	198	! get column / layer conversion factor
[1308]	199	DO k = 1, nlayer
[787]	200	DO i = 1, ngrid
[253]	201	l2c(i,k)=(pplev(i,k)-pplev(i,k+1))/g
[135]	202	ENDDO
	203	ENDDO
	204
	205	! Vertical loop (from top to bottom)
	206	! We carry the rain with us and calculate that added by warm/cold precipitation
	207	! processes and that subtracted by evaporation at each level.
[1308]	208	DO k = nlayer, 1, -1
[135]	209
	210	IF (evap_prec) THEN ! note no rneb dependence!
[787]	211	DO i = 1, ngrid
[135]	212	IF (zrfl(i) .GT.0.) THEN
[253]	213
[728]	214	if(zt(i,k).gt.Tsat(i,k))then
[863]	215	!! treat the case where all liquid water should boil
	216	zqev=MIN((zt(i,k)-Tsat(i,k))RCPDl2c(i,k)/RLVTT/ptimestep,zrfl(i))
[728]	217	zrfl(i)=MAX(zrfl(i)-zqev,0.)
[863]	218	d_q(i,k)=zqev/l2c(i,k)*ptimestep
[728]	219	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD
	220	else
[731]	221	zqev = MAX (0.0, (zqs(i,k)-q(i,k)))*l2c(i,k)/ptimestep !there was a bug here
[728]	222	zqevt= 2.0e-5*(1.0-q(i,k)/zqs(i,k)) & !default was 2.e-5
	223	sqrt(zrfl(i))l2c(i,k)/pplay(i,k)zt(i,k)R ! BC modif here
	224	zqevt = MAX (zqevt, 0.0)
	225	zqev = MIN (zqev, zqevt)
	226	zqev = MAX (zqev, 0.0)
	227	zrfln(i)= zrfl(i) - zqev
	228	zrfln(i)= max(zrfln(i),0.0)
[253]	229
[728]	230	d_q(i,k) = - (zrfln(i)-zrfl(i))/l2c(i,k)*ptimestep
	231	!d_t(i,k) = d_q(i,k) * RLVTT/RCPD!/(1.0+RVTMP2*q(i,k)) ! double BC modif here
	232	d_t(i,k) = - d_q(i,k) * RLVTT/RCPD ! was bugged!
	233	zrfl(i) = zrfln(i)
	234	end if
	235
[135]	236
[1283]	237	ENDIF ! of IF (zrfl(i) .GT.0.)
[135]	238	ENDDO
[1283]	239	ENDIF ! of IF (evap_prec)
[135]	240
[1283]	241	zoliq(1:ngrid) = 0.0
[135]	242
	243
[728]	244	if(precip_scheme.eq.1)then
[135]	245
[787]	246	DO i = 1, ngrid
[253]	247	ttemp = zt(i,k)
[650]	248	IF (ttemp .ge. T_h2O_ice_liq) THEN
[253]	249	lconvert=rainthreshold
	250	ELSEIF (ttemp .gt. t_crit) THEN
	251	lconvert=rainthreshold*(1.- t_crit/ttemp)
	252	lconvert=MAX(0.0,lconvert)
	253	ELSE
	254	lconvert=0.
	255	ENDIF
[135]	256
[253]	257
	258	IF (ql(i,k).gt.1.e-9) then
	259	zneb(i) = MAX(rneb(i,k), seuil_neb)
	260	IF ((ql(i,k)/zneb(i)).gt.lconvert)THEN ! precipitate!
[622]	261	d_ql(i,k) = -MAX((ql(i,k)-lconvert*zneb(i)),0.0)
[253]	262	zrfl(i) = zrfl(i) - d_ql(i,k)*l2c(i,k)/ptimestep
	263	ENDIF
	264	ENDIF
[135]	265	ENDDO
	266
[728]	267	elseif (precip_scheme.ge.2) then
	268
[787]	269	DO i = 1, ngrid
[135]	270	IF (rneb(i,k).GT.0.0) THEN
	271	zoliq(i) = ql(i,k)
[253]	272	zrho(i) = pplay(i,k) / ( zt(i,k) * R )
[135]	273	zdz(i) = (pplev(i,k)-pplev(i,k+1)) / (zrho(i)*g)
[650]	274	zfrac(i) = (zt(i,k)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
[135]	275	zfrac(i) = MAX(zfrac(i), 0.0)
	276	zfrac(i) = MIN(zfrac(i), 1.0)
	277	zneb(i) = MAX(rneb(i,k), seuil_neb)
	278	ENDIF
[731]	279	ENDDO
[135]	280
[731]	281	!recalculate liquid water particle radii
[1308]	282	call h2o_cloudrad(ngrid,nlayer,ql,reffh2oliq,reffh2oice)
[731]	283
[728]	284	SELECT CASE(precip_scheme)
	285	!precip scheme from Sundquist 78
	286	CASE(2)
	287
[135]	288	DO n = 1, ninter
[787]	289	DO i = 1, ngrid
[135]	290	IF (rneb(i,k).GT.0.0) THEN
[728]	291	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
[253]	292
[728]	293	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescale)) zoliq(i) &
	294	* (1.0-EXP(-(zoliq(i)/zneb(i)/cloud_sat)*2)) zfrac(i)
[135]	295	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	296	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
[731]	297	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
[135]	298	ztot(i) = zchau(i) + zfroi(i)
	299
	300	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	301	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	302	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
[253]	303
[135]	304	ENDIF
	305	ENDDO
	306	ENDDO
	307
[728]	308	!precip scheme modified from Sundquist 78 (in q**3)
	309	CASE(3)
	310
	311	DO n = 1, ninter
[787]	312	DO i = 1, ngrid
[728]	313	IF (rneb(i,k).GT.0.0) THEN
	314	! this is the ONLY place where zneb, precip_timescale and cloud_sat are used
	315
	316	zchau(i) = (ptimestep/(FLOAT(ninter)precip_timescalecloud_sat*2)) (zoliq(i)/zneb(i))**3
	317	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	318	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
[731]	319	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
[728]	320	ztot(i) = zchau(i) + zfroi(i)
	321
	322	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	323	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	324	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
	325
	326	ENDIF
	327	ENDDO
	328	ENDDO
	329
	330	!precip scheme modified from Boucher 95
	331	CASE(4)
	332
	333	DO n = 1, ninter
[787]	334	DO i = 1, ngrid
[728]	335	IF (rneb(i,k).GT.0.0) THEN
	336	! this is the ONLY place where zneb and c1 are used
	337
	338	zchau(i) = ptimestep/FLOAT(ninter) c1 zrho(i) &
	339	(zoliq(i)/zneb(i))2reffh2oliq(i,k)zneb(i) zfrac(i)
	340	zrhol(i) = zrho(i) * zoliq(i) / zneb(i)
	341	zfroi(i) = ptimestep/FLOAT(ninter)/zdz(i)*zoliq(i) &
[731]	342	fall2v(reffh2oice(i,k)) (1.0-zfrac(i)) ! zfroi behaves oddly...
[728]	343	ztot(i) = zchau(i) + zfroi(i)
	344
	345	IF (zneb(i).EQ.seuil_neb) ztot(i) = 0.0
	346	ztot(i) = MIN(MAX(ztot(i),0.0),zoliq(i))
	347	zoliq(i) = MAX(zoliq(i)-ztot(i), 0.0)
	348
	349	ENDIF
	350	ENDDO
	351	ENDDO
	352
	353	END SELECT ! precip_scheme
	354
[135]	355	! Change in cloud density and surface H2O values
[787]	356	DO i = 1, ngrid
[135]	357	IF (rneb(i,k).GT.0.0) THEN
[253]	358	d_ql(i,k) = (zoliq(i) - ql(i,k))!/ptimestep
	359	zrfl(i) = zrfl(i)+ MAX(ql(i,k)-zoliq(i),0.0)*l2c(i,k)/ptimestep
[135]	360	ENDIF
	361	ENDDO
	362
	363
[728]	364	endif ! if precip_scheme=1
	365
[1308]	366	ENDDO ! of DO k = nlayer, 1, -1
[135]	367
	368	! Rain or snow on the ground
[787]	369	DO i = 1, ngrid
[253]	370	if(zrfl(i).lt.0.0)then
	371	print*,'Droplets of negative rain are falling...'
	372	call abort
	373	endif
[650]	374	IF (t(i,1) .LT. T_h2O_ice_liq) THEN
[135]	375	dqssnow(i) = zrfl(i)
[253]	376	dqsrain(i) = 0.0
[135]	377	ELSE
[253]	378	dqssnow(i) = 0.0
[135]	379	dqsrain(i) = zrfl(i) ! liquid water = ice for now
	380	ENDIF
	381	ENDDO
	382
	383	! now subroutine -----> GCM variables
[1283]	384	if (evap_prec) then
[1308]	385	dqrain(1:ngrid,1:nlayer,i_vap)=d_q(1:ngrid,1:nlayer)/ptimestep
	386	d_t(1:ngrid,1:nlayer)=d_t(1:ngrid,1:nlayer)/ptimestep
[1283]	387	else
[1308]	388	dqrain(1:ngrid,1:nlayer,i_vap)=0.0
	389	d_t(1:ngrid,1:nlayer)=0.0
[1283]	390	endif
[1308]	391	dqrain(1:ngrid,1:nlayer,i_ice) = d_ql(1:ngrid,1:nlayer)/ptimestep
[135]	392
	393	end subroutine rain

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

Download in other formats:

Original Format