Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

verttransform.f90 @ 8

Last change on this file since 8 was 8, checked in by hasod, 11 years ago
Added parallel version of Flexpart91
File size: 20.3 KB

Line
1	!**********************************************************************
2	! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
3	! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
4	! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
5	! *
6	! This file is part of FLEXPART. *
7	! *
8	! FLEXPART is free software: you can redistribute it and/or modify *
9	! it under the terms of the GNU General Public License as published by*
10	! the Free Software Foundation, either version 3 of the License, or *
11	! (at your option) any later version. *
12	! *
13	! FLEXPART is distributed in the hope that it will be useful, *
14	! but WITHOUT ANY WARRANTY; without even the implied warranty of *
15	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
16	! GNU General Public License for more details. *
17	! *
18	! You should have received a copy of the GNU General Public License *
19	! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
20	!**********************************************************************
21
22	subroutine verttransform(n,uuh,vvh,wwh,pvh)
23	! i i i i i
24	!*****************************************************************************
25	! *
26	! This subroutine transforms temperature, dew point temperature and *
27	! wind components from eta to meter coordinates. *
28	! The vertical wind component is transformed from Pa/s to m/s using *
29	! the conversion factor pinmconv. *
30	! In addition, this routine calculates vertical density gradients *
31	! needed for the parameterization of the turbulent velocities. *
32	! *
33	! Author: A. Stohl, G. Wotawa *
34	! *
35	! 12 August 1996 *
36	! Update: 16 January 1998 *
37	! *
38	! Major update: 17 February 1999 *
39	! by G. Wotawa *
40	! *
41	! - Vertical levels for u, v and w are put together *
42	! - Slope correction for vertical velocity: Modification of calculation *
43	! procedure *
44	! *
45	!*****************************************************************************
46	! Changes, Bernd C. Krueger, Feb. 2001:
47	! Variables tth and qvh (on eta coordinates) from common block
48	!*****************************************************************************
49	! Sabine Eckhardt, March 2007
50	! added the variable cloud for use with scavenging - descr. in com_mod
51	!*****************************************************************************
52	! *
53	! Variables: *
54	! nx,ny,nz field dimensions in x,y and z direction *
55	! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition *
56	! uu(0:nxmax,0:nymax,nzmax,2) wind components in x-direction [m/s] *
57	! vv(0:nxmax,0:nymax,nzmax,2) wind components in y-direction [m/s] *
58	! ww(0:nxmax,0:nymax,nzmax,2) wind components in z-direction [deltaeta/s]*
59	! tt(0:nxmax,0:nymax,nzmax,2) temperature [K] *
60	! pv(0:nxmax,0:nymax,nzmax,2) potential voriticity (pvu) *
61	! ps(0:nxmax,0:nymax,2) surface pressure [Pa] *
62	! *
63	!*****************************************************************************
64
65	use par_mod
66	use com_mod
67	use cmapf_mod, only: cc2gll
68
69	implicit none
70
71	integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
72	integer :: rain_cloud_above,kz_inv
73	real :: f_qvsat,pressure
74	real :: rh,lsp,convp
75	real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
76	real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
77	real :: xlon,ylat,xlonr,dzdx,dzdy
78	real :: dzdx1,dzdx2,dzdy1,dzdy2
79	real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
80	real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
81	real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
82	real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
83	real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
84	real :: wzlev(nwzmax)
85	real,allocatable,dimension(:,:,:) :: uvwzlev
86	real,parameter :: const=r_air/ga
87
88	integer :: stat
89	logical :: init = .true.
90
91	allocate(uvwzlev(0:nxmax-1,0:nymax-1,nzmax),stat=stat)
92
93	!*************************************************************************
94	! If verttransform is called the first time, initialize heights of the *
95	! z levels in meter. The heights are the heights of model levels, where *
96	! u,v,T and qv are given, and of the interfaces, where w is given. So, *
97	! the vertical resolution in the z system is doubled. As reference point,*
98	! the lower left corner of the grid is used. *
99	! Unlike in the eta system, no difference between heights for u,v and *
100	! heights for w exists. *
101	!*************************************************************************
102
103
104	! do 897 kz=1,nuvz
105	! write (,) 'akz: ',akz(kz),'bkz',bkz(kz)
106	!897 continue
107
108	if (init) then
109
110	! Search for a point with high surface pressure (i.e. not above significant topography)
111	! Then, use this point to construct a reference z profile, to be used at all times
112	!*****************************************************************************
113
114	do jy=0,nymin1
115	do ix=0,nxmin1
116	if (ps(ix,jy,1,n).gt.100000.) then
117	ixm=ix
118	jym=jy
119	goto 3
120	endif
121	end do
122	end do
123	3 continue
124
125
126	tvold=tt2(ixm,jym,1,n)(1.+0.378ew(td2(ixm,jym,1,n))/ &
127	ps(ixm,jym,1,n))
128	pold=ps(ixm,jym,1,n)
129	height(1)=0.
130
131	do kz=2,nuvz
132	pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
133	tv=tth(ixm,jym,kz,n)(1.+0.608qvh(ixm,jym,kz,n))
134
135
136	! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
137	! upon the transformation to z levels. In order to save computer memory, this is
138	! not done anymore in the standard version. However, this option can still be
139	! switched on by replacing the following lines with those below, that are
140	! currently commented out.
141	! Note that two more changes are necessary in this subroutine below.
142	! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
143	!*****************************************************************************
144
145	if (abs(tv-tvold).gt.0.2) then
146	height(kz)= &
147	height(kz-1)+constlog(pold/pint) &
148	(tv-tvold)/log(tv/tvold)
149	else
150	height(kz)=height(kz-1)+ &
151	constlog(pold/pint)tv
152	endif
153
154	! Switch on following lines to use doubled vertical resolution
155	!*************************************************************
156	! if (abs(tv-tvold).gt.0.2) then
157	! height((kz-1)*2)=
158	! + height(max((kz-2)2,1))+constlog(pold/pint)*
159	! + (tv-tvold)/log(tv/tvold)
160	! else
161	! height((kz-1)2)=height(max((kz-2)2,1))+
162	! + constlog(pold/pint)tv
163	! endif
164	! End doubled vertical resolution
165
166	tvold=tv
167	pold=pint
168	end do
169
170
171	! Switch on following lines to use doubled vertical resolution
172	!*************************************************************
173	! do 7 kz=3,nz-1,2
174	! height(kz)=0.5*(height(kz-1)+height(kz+1))
175	! height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
176	! End doubled vertical resolution
177
178
179	! Determine highest levels that can be within PBL
180	!************************************************
181
182	do kz=1,nz
183	if (height(kz).gt.hmixmax) then
184	nmixz=kz
185	goto 9
186	endif
187	end do
188	9 continue
189
190	! Do not repeat initialization of the Cartesian z grid
191	!*****************************************************
192
193	init=.false.
194
195	endif
196
197
198	! Loop over the whole grid
199	!*************************
200
201	do jy=0,nymin1
202	do ix=0,nxmin1
203	tvold=tt2(ix,jy,1,n)(1.+0.378ew(td2(ix,jy,1,n))/ &
204	ps(ix,jy,1,n))
205	pold=ps(ix,jy,1,n)
206	uvzlev(1)=0.
207	wzlev(1)=0.
208	rhoh(1)=pold/(r_air*tvold)
209
210
211	! Compute heights of eta levels
212	!******************************
213
214	do kz=2,nuvz
215	pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
216	tv=tth(ix,jy,kz,n)(1.+0.608qvh(ix,jy,kz,n))
217	rhoh(kz)=pint/(r_air*tv)
218
219	if (abs(tv-tvold).gt.0.2) then
220	uvzlev(kz)=uvzlev(kz-1)+constlog(pold/pint) &
221	(tv-tvold)/log(tv/tvold)
222	else
223	uvzlev(kz)=uvzlev(kz-1)+constlog(pold/pint)tv
224	endif
225
226	tvold=tv
227	pold=pint
228	end do
229
230
231	do kz=2,nwz-1
232	wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
233	end do
234	wzlev(nwz)=wzlev(nwz-1)+ &
235	uvzlev(nuvz)-uvzlev(nuvz-1)
236
237	uvwzlev(ix,jy,1)=0.0
238	do kz=2,nuvz
239	uvwzlev(ix,jy,kz)=uvzlev(kz)
240	end do
241
242	! Switch on following lines to use doubled vertical resolution
243	! Switch off the three lines above.
244	!*************************************************************
245	!22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
246	! do 23 kz=2,nwz
247	!23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
248	! End doubled vertical resolution
249
250	! pinmconv=(h2-h1)/(p2-p1)
251
252	pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
253	((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
254	(aknew(1)+bknew(1)*ps(ix,jy,1,n)))
255	do kz=2,nz-1
256	pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
257	((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
258	(aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
259	end do
260	pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
261	((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
262	(aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
263
264	! Levels, where u,v,t and q are given
265	!************************************
266
267	uu(ix,jy,1,n)=uuh(ix,jy,1)
268	vv(ix,jy,1,n)=vvh(ix,jy,1)
269	tt(ix,jy,1,n)=tth(ix,jy,1,n)
270	qv(ix,jy,1,n)=qvh(ix,jy,1,n)
271	pv(ix,jy,1,n)=pvh(ix,jy,1)
272	rho(ix,jy,1,n)=rhoh(1)
273	uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
274	vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
275	tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
276	qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
277	pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
278	rho(ix,jy,nz,n)=rhoh(nuvz)
279	kmin=2
280	do iz=2,nz-1
281	do kz=kmin,nuvz
282	if(height(iz).gt.uvzlev(nuvz)) then
283	uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
284	vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
285	tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
286	qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
287	pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
288	rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
289	goto 30
290	endif
291	if ((height(iz).gt.uvzlev(kz-1)).and. &
292	(height(iz).le.uvzlev(kz))) then
293	dz1=height(iz)-uvzlev(kz-1)
294	dz2=uvzlev(kz)-height(iz)
295	dz=dz1+dz2
296	uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)dz2+uuh(ix,jy,kz)dz1)/dz
297	vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)dz2+vvh(ix,jy,kz)dz1)/dz
298	tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
299	+tth(ix,jy,kz,n)*dz1)/dz
300	qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
301	+qvh(ix,jy,kz,n)*dz1)/dz
302	pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)dz2+pvh(ix,jy,kz)dz1)/dz
303	rho(ix,jy,iz,n)=(rhoh(kz-1)dz2+rhoh(kz)dz1)/dz
304	kmin=kz
305	goto 30
306	endif
307	end do
308	30 continue
309	end do
310
311
312	! Levels, where w is given
313	!*************************
314
315	ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
316	ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
317	kmin=2
318	do iz=2,nz
319	do kz=kmin,nwz
320	if ((height(iz).gt.wzlev(kz-1)).and. &
321	(height(iz).le.wzlev(kz))) then
322	dz1=height(iz)-wzlev(kz-1)
323	dz2=wzlev(kz)-height(iz)
324	dz=dz1+dz2
325	ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)pinmconv(kz-1)dz2 &
326	+wwh(ix,jy,kz)pinmconv(kz)dz1)/dz
327	kmin=kz
328	goto 40
329	endif
330	end do
331	40 continue
332	end do
333
334	! Compute density gradients at intermediate levels
335	!*************************************************
336
337	drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
338	(height(2)-height(1))
339	do kz=2,nz-1
340	drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
341	(height(kz+1)-height(kz-1))
342	end do
343	drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
344
345	end do
346	end do
347
348
349	!****************************************************************
350	! Compute slope of eta levels in windward direction and resulting
351	! vertical wind correction
352	!****************************************************************
353
354	do jy=1,ny-2
355	do ix=1,nx-2
356
357	kmin=2
358	do iz=2,nz-1
359
360	ui=uu(ix,jy,iz,n)dxconst/cos((real(jy)dy+ylat0)*pi180)
361	vi=vv(ix,jy,iz,n)*dyconst
362
363	do kz=kmin,nz
364	if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
365	(height(iz).le.uvwzlev(ix,jy,kz))) then
366	dz1=height(iz)-uvwzlev(ix,jy,kz-1)
367	dz2=uvwzlev(ix,jy,kz)-height(iz)
368	dz=dz1+dz2
369	kl=kz-1
370	klp=kz
371	kmin=kz
372	goto 47
373	endif
374	end do
375
376	47 ix1=ix-1
377	jy1=jy-1
378	ixp=ix+1
379	jyp=jy+1
380
381	dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
382	dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
383	dzdx=(dzdx1dz2+dzdx2dz1)/dz
384
385	dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
386	dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
387	dzdy=(dzdy1dz2+dzdy2dz1)/dz
388
389	ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdxui+dzdyvi)
390
391	end do
392
393	end do
394	end do
395
396
397	! If north pole is in the domain, calculate wind velocities in polar
398	! stereographic coordinates
399	!*******************************************************************
400
401	if (nglobal) then
402	do jy=int(switchnorthg)-2,nymin1
403	ylat=ylat0+real(jy)*dy
404	do ix=0,nxmin1
405	xlon=xlon0+real(ix)*dx
406	do iz=1,nz
407	call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
408	vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
409	vvpol(ix,jy,iz,n))
410	end do
411	end do
412	end do
413
414
415	do iz=1,nz
416
417	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
418	!
419	! AMSnauffer Nov 18 2004 Added check for case vv=0
420	!
421	xlon=xlon0+real(nx/2-1)*dx
422	xlonr=xlon*pi/180.
423	ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
424	vv(nx/2-1,nymin1,iz,n)**2)
425	if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
426	ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
427	vv(nx/2-1,nymin1,iz,n))-xlonr
428	else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
429	ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
430	vv(nx/2-1,nymin1,iz,n))-xlonr
431	else
432	ddpol=pi/2-xlonr
433	endif
434	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
435	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
436
437	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
438	xlon=180.0
439	xlonr=xlon*pi/180.
440	ylat=90.0
441	uuaux=-ffpol*sin(xlonr+ddpol)
442	vvaux=-ffpol*cos(xlonr+ddpol)
443	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
444	vvpolaux)
445
446	jy=nymin1
447	do ix=0,nxmin1
448	uupol(ix,jy,iz,n)=uupolaux
449	vvpol(ix,jy,iz,n)=vvpolaux
450	end do
451	end do
452
453
454	! Fix: Set W at pole to the zonally averaged W of the next equator-
455	! ward parallel of latitude
456
457	do iz=1,nz
458	wdummy=0.
459	jy=ny-2
460	do ix=0,nxmin1
461	wdummy=wdummy+ww(ix,jy,iz,n)
462	end do
463	wdummy=wdummy/real(nx)
464	jy=nymin1
465	do ix=0,nxmin1
466	ww(ix,jy,iz,n)=wdummy
467	end do
468	end do
469
470	endif
471
472
473	! If south pole is in the domain, calculate wind velocities in polar
474	! stereographic coordinates
475	!*******************************************************************
476
477	if (sglobal) then
478	do jy=0,int(switchsouthg)+3
479	ylat=ylat0+real(jy)*dy
480	do ix=0,nxmin1
481	xlon=xlon0+real(ix)*dx
482	do iz=1,nz
483	call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
484	vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
485	vvpol(ix,jy,iz,n))
486	end do
487	end do
488	end do
489
490	do iz=1,nz
491
492	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
493	!
494	! AMSnauffer Nov 18 2004 Added check for case vv=0
495	!
496	xlon=xlon0+real(nx/2-1)*dx
497	xlonr=xlon*pi/180.
498	ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
499	vv(nx/2-1,0,iz,n)**2)
500	if (vv(nx/2-1,0,iz,n).lt.0.) then
501	ddpol=atan(uu(nx/2-1,0,iz,n)/ &
502	vv(nx/2-1,0,iz,n))+xlonr
503	else if (vv(nx/2-1,0,iz,n).gt.0.) then
504	ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
505	vv(nx/2-1,0,iz,n))+xlonr
506	else
507	ddpol=pi/2-xlonr
508	endif
509	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
510	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
511
512	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
513	xlon=180.0
514	xlonr=xlon*pi/180.
515	ylat=-90.0
516	uuaux=+ffpol*sin(xlonr-ddpol)
517	vvaux=-ffpol*cos(xlonr-ddpol)
518	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
519	vvpolaux)
520
521	jy=0
522	do ix=0,nxmin1
523	uupol(ix,jy,iz,n)=uupolaux
524	vvpol(ix,jy,iz,n)=vvpolaux
525	end do
526	end do
527
528
529	! Fix: Set W at pole to the zonally averaged W of the next equator-
530	! ward parallel of latitude
531
532	do iz=1,nz
533	wdummy=0.
534	jy=1
535	do ix=0,nxmin1
536	wdummy=wdummy+ww(ix,jy,iz,n)
537	end do
538	wdummy=wdummy/real(nx)
539	jy=0
540	do ix=0,nxmin1
541	ww(ix,jy,iz,n)=wdummy
542	end do
543	end do
544	endif
545
546
547	!write (,) 'initializing clouds, n:',n,nymin1,nxmin1,nz
548	! create a cloud and rainout/washout field, clouds occur where rh>80%
549	! total cloudheight is stored at level 0
550	do jy=0,nymin1
551	do ix=0,nxmin1
552	rain_cloud_above=0
553	lsp=lsprec(ix,jy,1,n)
554	convp=convprec(ix,jy,1,n)
555	cloudsh(ix,jy,n)=0
556	do kz_inv=1,nz-1
557	kz=nz-kz_inv+1
558	pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
559	rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
560	clouds(ix,jy,kz,n)=0
561	if (rh.gt.0.8) then ! in cloud
562	if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
563	rain_cloud_above=1
564	cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
565	int(height(kz)-height(kz-1))
566	if (lsp.ge.convp) then
567	clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
568	else
569	clouds(ix,jy,kz,n)=2 ! convp dominated rainout
570	endif
571	else ! no precipitation
572	clouds(ix,jy,kz,n)=1 ! cloud
573	endif
574	else ! no cloud
575	if (rain_cloud_above.eq.1) then ! scavenging
576	if (lsp.ge.convp) then
577	clouds(ix,jy,kz,n)=5 ! lsp dominated washout
578	else
579	clouds(ix,jy,kz,n)=4 ! convp dominated washout
580	endif
581	endif
582	endif
583	end do
584	end do
585	end do
586
587	!do 102 kz=1,nuvz
588	!write(an,'(i02)') kz+10
589	!write(,) nuvz,nymin1,nxmin1,'--',an,'--'
590	!open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
591	!do 101 jy=0,nymin1
592	! write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
593	!101 continue
594	! close(4)
595	!102 continue
596
597	! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
598	! do 103 jy=0,nymin1
599	! write (4,*)
600	!+ (height(kz),kz=1,nuvz)
601	!103 continue
602	! close(4)
603
604	!open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
605	! do 104 jy=0,nymin1
606	! write (4,*)
607	!+ (r_airtt(ix,jy,1,n)rho(ix,jy,1,n),ix=1,nxmin1)
608	!104 continue
609	! close(4)
610
611	deallocate(uvwzlev)
612
613	end subroutine verttransform

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

Download in other formats:

Original Format