Context Navigation

moistadj.F90 @ 735

Last change on this file since 735 was 728, checked in by jleconte, 13 years ago

18/07/2012 == JL

New water cycle scheme:
- largescale now in F90. Robustness increased by i) including evap inside largescale ii) computing the

condensed water amount iteratively

same improvements in moistadj.
Water thermodynamical data and saturation curves centralized in module watercommn_h
- The saturation curves used are now Tetens formula as they are analyticaly inversible (Ts(P)-> Ps(T)).

New saturation curve yields very good agreement with the former one.

Saturation curves are now generalized for arbitrary water amount (not just q<<1)
The old watersat should be removed soon.

The effect of water vapor on total (surface) pressure can be taken into account by setting

mass_redistrib=.true. in callphys.def (routine mass_redistribution inspired from co2_condense in martian
model but with a different scheme as many routines evaporate/condense water vapor).

New cloud and precipitation scheme (JL + BC):
- The default recovery assumption for computing the total cloud fraction has been changed (total random gave too

large cloud fractions). See totalcloudfrac.F90 for details and to change this.

Totalcloudfraction now set the total cloud fraction to the fraction of the

optically thickest cloud and totalcloudfrac is thus called in aeropacity.

Only the total cloud fraction is used to compute optical depth in aeropacity (no more effective

optical depth with exponential formula).

4 precipitation schemes are now available (see rain.F90 for details). The choice can be made using precip_scheme

in callphys.def. Usage of the more physically based model of Boucher et al 95 (precip_scheme=4) is recommended.
default behavior is set to the former "simple scheme" (precip_scheme=1).

See rain.f90 to determine the parameter to be defined in callphys.def as a function of the precipitation scheme used.

Physiq.F90 now written in a matricial (more F90) way.
Radii (H2O and CO2 cloud particles, aerosols, duts, ...) calculations now centralized in module radii_mod.F90

and work with the new aerosol scheme implemented by Laura K. Some inconsistency may remain in callsedim.

Implementation compiled with ifort and pgf90.
gcm.e runs in Earth and Early Mars case with CO2 and H2O cycle + dust.

File size: 10.2 KB

Rev	Line
[728]	1	subroutine moistadj(pt, pq, pdq, pplev, pplay, pdtmana, pdqmana, ptimestep, rneb)
[135]	2
[728]	3	use watercommon_h, only: T_h2O_ice_liq, RLVTT, RCPD, Psat_water, Lcpdqsat_water
[135]	4
	5	implicit none
	6
	7
	8	!=====================================================================
	9	!
	10	! Purpose
	11	! -------
	12	! Calculates moist convective adjustment by the method of Manabe.
	13	!
	14	! Authors
	15	! -------
	16	! Adapted from the LMDTERRE code by R. Wordsworth (2010)
	17	! Original author Z. X. Li (1993)
	18	!
	19	!=====================================================================
	20
	21	#include "dimensions.h"
	22	#include "dimphys.h"
	23	#include "tracer.h"
	24	#include "comcstfi.h"
	25
[728]	26	REAL pt(ngridmx,nlayermx) ! temperature (K)
	27	REAL pq(ngridmx,nlayermx,nqmx) ! tracer (kg/kg)
[135]	28	REAL pdq(ngridmx,nlayermx,nqmx)
	29
[728]	30	REAL pdqmana(ngridmx,nlayermx,nqmx) ! tendency of tracers (kg/kg.s-1)
	31	REAL pdtmana(ngridmx,nlayermx) ! temperature increment
[135]	32
[728]	33	! local variables
	34	REAL zt(ngridmx,nlayermx) ! temperature (K)
	35	REAL zq(ngridmx,nlayermx) ! humidite specifique (kg/kg)
[135]	36	REAL pplev(ngridmx,nlayermx+1) ! pression a inter-couche (Pa)
	37	REAL pplay(ngridmx,nlayermx) ! pression au milieu de couche (Pa)
	38
	39	REAL d_t(ngridmx,nlayermx) ! temperature increment
	40	REAL d_q(ngridmx,nlayermx) ! incrementation pour vapeur d'eau
	41	REAL d_ql(ngridmx,nlayermx) ! incrementation pour l'eau liquide
	42	REAL rneb(ngridmx,nlayermx) ! cloud fraction
	43	REAL ptimestep
	44
	45	! REAL t_coup
	46	! PARAMETER (t_coup=234.0)
	47	REAL seuil_vap
	48	PARAMETER (seuil_vap=1.0E-10)
	49
	50	! Local variables
[728]	51	INTEGER i, k, iq, nn
	52	INTEGER, PARAMETER :: niter = 6
[135]	53	INTEGER k1, k1p, k2, k2p
	54	LOGICAL itest(ngridmx)
	55	REAL delta_q(ngridmx, nlayermx)
	56	REAL cp_new_t(nlayermx)
	57	REAL cp_delta_t(nlayermx)
	58	REAL v_cptj(nlayermx), v_cptjk1, v_ssig
[728]	59	REAL v_cptt(ngridmx,nlayermx), v_p, v_t, zpsat
	60	REAL zqs(ngridmx,nlayermx), zdqs(ngridmx,nlayermx)
[135]	61	REAL zq1(ngridmx), zq2(ngridmx)
	62	REAL gamcpdz(ngridmx,2:nlayermx)
	63	REAL zdp, zdpm
	64
	65	REAL zsat ! super-saturation
	66	REAL zflo ! flotabilite
	67
	68	REAL local_q(ngridmx,nlayermx),local_t(ngridmx,nlayermx)
	69
	70	REAL zdelta, zcor, zcvm5
	71
	72	REAL dEtot, dqtot, masse ! conservation diagnostics
	73	real dL1tot, dL2tot
	74
	75	! Indices of water vapour and water ice tracers
	76	INTEGER,SAVE :: i_h2o=0 ! water vapour
	77	INTEGER,SAVE :: i_ice=0 ! water ice
	78
	79	LOGICAL firstcall
	80	SAVE firstcall
	81
	82	DATA firstcall /.TRUE./
	83
	84	IF (firstcall) THEN
	85
	86	i_h2o=igcm_h2o_vap
	87	i_ice=igcm_h2o_ice
	88
	89	write(,) "rain: i_ice=",i_ice
	90	write(,) " i_h2o=",i_h2o
	91
	92	firstcall = .FALSE.
	93	ENDIF
	94
	95	! GCM -----> subroutine variables
[728]	96	zq(1:ngridmx,1:nlayermx) = pq(1:ngridmx,1:nlayermx,i_h2o)+ pdq(1:ngridmx,1:nlayermx,i_h2o)*ptimestep
	97	zt(1:ngridmx,1:nlayermx) = pt(1:ngridmx,1:nlayermx)
	98	pdqmana(1:ngridmx,1:nlayermx,1:nqmx)=0.0
	99
[135]	100	DO k = 1, nlayermx
[728]	101	DO i = 1, ngridmx
	102	if(zq(i,k).lt.0.)then
	103	zq(i,k)=0.0
[135]	104	endif
[728]	105	ENDDO
[135]	106	ENDDO
[728]	107
	108	local_q(1:ngridmx,1:nlayermx) = zq(1:ngridmx,1:nlayermx)
	109	local_t(1:ngridmx,1:nlayermx) = zt(1:ngridmx,1:nlayermx)
	110	rneb(1:ngridmx,1:nlayermx) = 0.0
	111	d_ql(1:ngridmx,1:nlayermx) = 0.0
	112	d_t(1:ngridmx,1:nlayermx) = 0.0
	113	d_q(1:ngridmx,1:nlayermx) = 0.0
[135]	114
	115	! Calculate v_cptt
	116	DO k = 1, nlayermx
	117	DO i = 1, ngridmx
	118	v_cptt(i,k) = RCPD * local_t(i,k)
[728]	119	v_t = MAX(local_t(i,k),15.)
[135]	120	v_p = pplay(i,k)
	121
[728]	122	call Psat_water(v_t,v_p,zpsat,zqs(i,k))
	123	call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
[135]	124	ENDDO
	125	ENDDO
	126
	127	! Calculate Gamma * Cp * dz: (gamma is the critical gradient)
	128	DO k = 2, nlayermx
	129	DO i = 1, ngridmx
	130	zdp = pplev(i,k)-pplev(i,k+1)
	131	zdpm = pplev(i,k-1)-pplev(i,k)
[728]	132	gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
	133	+ RLVTT /(zdpm+zdp) * (zqs(i,k-1)zdpm + zqs(i,k)zdp) ) &
	134	* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
	135	/ (1.0+ (zdqs(i,k-1)zdpm + zdqs(i,k)zdp)/(zdpm+zdp) )
[135]	136	ENDDO
	137	ENDDO
	138
	139	!------------------------------------ modification of unstable profile
	140	DO 9999 i = 1, ngridmx
	141	itest(i) = .FALSE.
	142
	143	! print*,'we in the loop'
	144	! stop
	145
	146	k1 = 0
	147	k2 = 1
	148
	149	810 CONTINUE ! look for k1, the base of the column
	150	k2 = k2 + 1
	151	IF (k2 .GT. nlayermx) GOTO 9999
	152	zflo = v_cptt(i,k2-1) - v_cptt(i,k2) - gamcpdz(i,k2)
[728]	153	zsat=(local_q(i,k2-1)-zqs(i,k2-1))*(pplev(i,k2-1)-pplev(i,k2)) &
	154	+(local_q(i,k2)-zqs(i,k2))*(pplev(i,k2)-pplev(i,k2+1))
[135]	155
	156	IF ( zflo.LE.0.0 .OR. zsat.LE.0.0 ) GOTO 810
	157	k1 = k2 - 1
	158	itest(i) = .TRUE.
	159
	160	820 CONTINUE !! look for k2, the top of the column
	161	IF (k2 .EQ. nlayermx) GOTO 821
	162	k2p = k2 + 1
[728]	163	zsat=zsat+(pplev(i,k2p)-pplev(i,k2p+1))*(local_q(i,k2p)-zqs(i,k2p))
[135]	164	zflo = v_cptt(i,k2p-1) - v_cptt(i,k2p) - gamcpdz(i,k2p)
	165
	166	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 821
	167	k2 = k2p
	168	GOTO 820
	169	821 CONTINUE
	170
	171	!------------------------------------------------------ local adjustment
	172	830 CONTINUE ! actual adjustment
[728]	173	Do nn=1,niter
[135]	174	v_cptj(k1) = 0.0
	175	zdp = pplev(i,k1)-pplev(i,k1+1)
[728]	176	v_cptjk1 = ( (1.0+zdqs(i,k1))(v_cptt(i,k1)+v_cptj(k1)) + RLVTT(local_q(i,k1)-zqs(i,k1)) ) * zdp
	177	v_ssig = zdp * (1.0+zdqs(i,k1))
[135]	178
	179	k1p = k1 + 1
	180	DO k = k1p, k2
	181	zdp = pplev(i,k)-pplev(i,k+1)
	182	v_cptj(k) = v_cptj(k-1) + gamcpdz(i,k)
[728]	183	v_cptjk1 = v_cptjk1 + zdp * ( (1.0+zdqs(i, k))(v_cptt(i,k)+v_cptj(k)) + RLVTT(local_q(i,k)-zqs(i,k)) )
	184	v_ssig = v_ssig + zdp *(1.0+zdqs(i,k))
[135]	185	ENDDO
	186
	187
	188	! this right here is where the adjustment is done???
	189	DO k = k1, k2
	190	cp_new_t(k) = v_cptjk1/v_ssig - v_cptj(k)
	191	cp_delta_t(k) = cp_new_t(k) - v_cptt(i,k)
[728]	192	v_cptt(i,k)=cp_new_t(k)
	193	local_q(i,k) = zqs(i,k) + zdqs(i,k)*cp_delta_t(k)/RLVTT
[135]	194	local_t(i,k) = cp_new_t(k) / RCPD
[253]	195
[728]	196	v_t = MAX(local_t(i,k),15.)
	197	v_p = pplay(i,k)
	198
	199	call Psat_water(v_t,v_p,zpsat,zqs(i,k))
	200	call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
	201
	202
	203
	204	! print*,'i,k,zqs,cpdT=',i,k,zqs(i,k),cp_delta_t(k)
[135]	205	ENDDO
[728]	206	Enddo ! nn=1,niter
[135]	207
[253]	208
[135]	209	!--------------------------------------------------- sounding downwards
	210	! -- we refine the prognostic variables in
	211	! -- the layer about to be adjusted
	212
[728]	213	! DO k = k1, k2
	214	! v_cptt(i,k) = RCPD * local_t(i,k)
	215	! v_t = local_t(i,k)
	216	! v_p = pplay(i,k)
	217	!
	218	! call Psat_water(v_t,v_p,zpsat,zqs(i,k))
	219	! call Lcpdqsat_water(v_t,v_p,zpsat,zqs(i,k),zdqs(i,k))
	220	! ENDDO
[135]	221
	222	DO k = 2, nlayermx
	223	zdpm = pplev(i,k-1) - pplev(i,k)
	224	zdp = pplev(i,k) - pplev(i,k+1)
[728]	225	gamcpdz(i,k) = ( ( R/RCPD /(zdpm+zdp) * (v_cptt(i,k-1)zdpm + v_cptt(i,k)zdp) &
	226	+ RLVTT /(zdpm+zdp) * (zqs(i,k-1)zdpm + zqs(i,k)zdp) ) &
	227	* (pplay(i,k-1)-pplay(i,k)) / pplev(i,k) ) &
	228	/ (1.0+ (zdqs(i,k-1)zdpm + zdqs(i,k)zdp)/(zdpm+zdp) )
[135]	229	ENDDO
	230
	231	! Test to see if we've reached the bottom
	232
	233	IF (k1 .EQ. 1) GOTO 841 ! yes we have!
	234	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
[728]	235	zsat=(local_q(i,k1-1)-zqs(i,k1-1))*(pplev(i,k1-1)-pplev(i,k1)) &
	236	+ (local_q(i,k1)-zqs(i,k1))*(pplev(i,k1)-pplev(i,k1+1))
[135]	237	IF (zflo.LE.0.0 .OR. zsat.LE.0.0) GOTO 841 ! yes we have!
	238
	239	840 CONTINUE
	240	k1 = k1 - 1
	241	IF (k1 .EQ. 1) GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
[728]	242	zsat = zsat + (local_q(i,k1-1)-zqs(i,k1-1)) &
[135]	243	*(pplev(i,k1-1)-pplev(i,k1))
	244	zflo = v_cptt(i,k1-1) - v_cptt(i,k1) - gamcpdz(i,k1)
	245	IF (zflo.GT.0.0 .AND. zsat.GT.0.0) THEN
	246	GOTO 840
	247	ELSE
	248	GOTO 830 ! GOTO 820 (a tester, Z.X.Li, mars 1995)
	249	ENDIF
	250	841 CONTINUE
	251
	252	GOTO 810 ! look for other layers higher up
	253
	254	9999 CONTINUE ! loop over all the points
	255
	256	! print*,'k1=',k1
	257	! print*,'k2=',k2
	258
	259	! print*,'local_t=',local_t
	260	! print*,'v_cptt=',v_cptt
	261	! print*,'gamcpdz=',gamcpdz
	262
	263	!-----------------------------------------------------------------------
	264	! Determine the cloud fraction (hypothese: la nebulosite a lieu
	265	! a l'endroit ou la vapeur d'eau est diminuee par l'ajustement):
	266
	267	DO k = 1, nlayermx
	268	DO i = 1, ngridmx
	269	IF (itest(i)) THEN
[728]	270	delta_q(i,k) = local_q(i,k) - zq(i,k)
[135]	271	IF (delta_q(i,k).LT.0.) rneb(i,k) = 1.0
	272	ENDIF
	273	ENDDO
	274	ENDDO
	275
	276	! Distribuer l'eau condensee en eau liquide nuageuse (hypothese:
	277	! l'eau liquide est distribuee aux endroits ou la vapeur d'eau
	278	! diminue et d'une maniere proportionnelle a cet diminution):
	279
	280	DO i = 1, ngridmx
	281	IF (itest(i)) THEN
	282	zq1(i) = 0.0
	283	zq2(i) = 0.0
	284	ENDIF
	285	ENDDO
	286	DO k = 1, nlayermx
	287	DO i = 1, ngridmx
	288	IF (itest(i)) THEN
	289	zdp = pplev(i,k)-pplev(i,k+1)
	290	zq1(i) = zq1(i) - delta_q(i,k) * zdp
	291	zq2(i) = zq2(i) - MIN(0.0, delta_q(i,k)) * zdp
	292	ENDIF
	293	ENDDO
	294	ENDDO
	295	DO k = 1, nlayermx
	296	DO i = 1, ngridmx
	297	IF (itest(i)) THEN
[728]	298	IF (zq2(i).NE.0.0) d_ql(i,k) = - MIN(0.0,delta_q(i,k))*zq1(i)/zq2(i)
[135]	299	ENDIF
	300	ENDDO
	301	ENDDO
	302
[253]	303	! print*,'local_q BEFORE=',local_q
	304
[135]	305	DO k = 1, nlayermx
	306	DO i = 1, ngridmx
	307	local_q(i, k) = MAX(local_q(i, k), seuil_vap)
	308	ENDDO
	309	ENDDO
	310
	311	DO k = 1, nlayermx
	312	DO i = 1, ngridmx
[728]	313	d_t(i,k) = local_t(i,k) - zt(i,k)
	314	d_q(i,k) = local_q(i,k) - zq(i,k)
[135]	315	ENDDO
	316	ENDDO
	317
	318	! now subroutine -----> GCM variables
	319	DO k = 1, nlayermx
	320	DO i = 1, ngridmx
	321
[728]	322	pdtmana(i,k) = d_t(i,k)/ptimestep
	323	pdqmana(i,k,i_h2o) = d_q(i,k)/ptimestep
	324	pdqmana(i,k,i_ice) = d_ql(i,k)/ptimestep
[135]	325
	326	ENDDO
	327	ENDDO
	328
[253]	329
[135]	330	RETURN
[253]	331	END

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

Download in other formats:

Original Format