Context Navigation

moistadj.F90 @ 848

Last change on this file since 848 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: 10.2 KB

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

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

Download in other formats:

Original Format