Context Navigation

rain.F90 @ 832

Last change on this file since 832 was 787, checked in by aslmd, 12 years ago

LMDZ.GENERIC. (Sorry for long text but this is a quite major commit)

Paving the path for parallel computations. And moving towards a more flexible code.

Automatic allocation is used within all routines in phystd. No further mention to ngridmx and nqmx.

ngridmx and nqmx are still used in LMDZ.GENERIC in the dyn3d part
if the LMDZ4/LMDZ5 dynamical core is used, there is no more fixed dimensions ngridmx and nqmx --> a fully flexible parallel implementation is now possible (e.g. no need to recompile when changing numbers of processors)

The important stuff :

Compilation checked with ifort. OK with and without debug mode. No errors. Checked for: gcm, newstart, rcm1d, kcm1d
RUN GCM: Running an Earth test case. Comparison with previous revision --> debug mode : perfect match. bit by bit (diff command). checked with plots --> O1 mode : close match (checked with plots) --> O2 mode : sometimes up to 0.5 K departure.... BUT in this new version O2 and O1 are quite close while in previous version O1 and O2 differed by about, well, typically 0.5 K (pictures available on request)
RUN NEWSTART : perfect match (bit-by-bit) in either debug or normal mode.
RUN RCM1D : perfect match in normal mode.
RUN KCM1D : not tested (I don't know what is the use of kcm1d)

List of main changes :

Additional arguments to some subroutines (ngrid and nq)
F77 include strategy is obsolete and replaced by F90 module strategy In this new strategy arrays are allocatable and allocated once at first use This has to be done for all common featuring arrays defined with ngridmx or nqmx

surfdat.h >> surfdat_h.F90
tracer.h >> tracer_h.F90
comsaison.h >> comsaison_h.F90
comgeomfi.h >> comgeomfi_h.F90
comsoil.h >> comsoil_h.F90
comdiurn.h >> comdiurn_h.F90
fisice.h >> DELETED. was not used. probably a fossil.
watercap.h >> DELETED. variable put in surfdat_h.F90

F77 'save' strategy is obsolete and replaced by F90 'allocatable save' strategy (see previous point and e.g. new version of physiq.F90)
Suppressing any mention to advtrac.h which is a common in the dynamics and needs nqmx This was easily solved by adding an argument with tracer names, coming from the dynamics This is probably not a definitive solution, ... but this allows for generic physics to work easily with either LMDZ.GENERIC or LMDZ dynamical cores
Removing consistency tests between nq and nqmx ; and ngrid and ngridmx. No use now!
Adaptation of rcm1d, kcm1d, newstart given above-mentioned changes

A note on phyetat0 and soil_setting:

Now written so that a slice of horizontal size 'ngrid' starting at grid point 'cursor' is read in startfi.nc 'cursor' is defined in dimphys.h and initialized by inifis (or in newstart) this is useful for parallel computations. default behavior is the usual one : sequential runs, cursor is 1, size ngrid is the whole global domain

A note on an additional change :

nueffrad is now an argument to callcorrk as is the case for reffrad both are saved in physiq this is for consistency and lisibility (previously nueffrad was saved in callcorrk) ... but there is a call to a function which modifies nueffrad in physiq ... previously this was not modifying nueffrad (although it was quite cumbersome to detect this) ... to be conservative I kept this behaviour and highlighted it with an array nueffrad_dummy ... I added a comment because someone might want to change this

File size: 13.7 KB

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

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

Download in other formats:

Original Format