Context Navigation

verttransform_ecmwf.f90 @ 8b3d324

univie

Last change on this file since 8b3d324 was 8b3d324, checked in by pesei <petra seibert at univie ac at>, 6 years ago

comment out print statements

comment out print statements for the reference height, adjust the i/o
comment line

Property mode set to 100644

File size: 32.6 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_ecmwf(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	! *
39	!*****************************************************************************
40	! CHANGES *
41	! Major update: 17 February 1999 *
42	! by G. Wotawa *
43	! *
44	! - Vertical levels for u, v and w are put together *
45	! - Slope correction for vertical velocity: Modification of calculation *
46	! procedure *
47	! *
48	! Bernd C. Krueger, Feb. 2001:
49	! Variables tth and qvh (on eta coordinates) from common block
50	! *
51	! Sabine Eckhardt, March 2007:
52	! added the variable cloud for use with scavenging - descr. in com_mod
53	! *
54	! Who? When? *
55	! Unified ECMWF and GFS builds
56	! Marian Harustak, 12.5.2017
57	! - Renamed from verttransform to verttransform_ecmwf
58	! *
59	! Don Morton, 2017-05-30:
60	! modification of a bug in ew. Don Morton (CTBTO project) *
61	! *
62	! undocumented modifications by NILU for v10 *
63	! *
64	! Petra Seibert, 2018-06-13: *
65	! - fix bug of ticket:140 (find reference position for vertical grid) *
66	! - put back SAVE attribute for INIT, just to be safe *
67	! - minor changes, most of them just cosmetics *
68	! for details see changelog.txt *
69	! *
70	! ****************************************************************************
71
72	!*****************************************************************************
73	! *
74	! Variables: *
75	! nx,ny,nz field dimensions in x,y and z direction *
76	! clouds(0:nxmax,0:nymax,0:nzmax,numwfmem) cloud field for wet deposition *
77	! uu(0:nxmax,0:nymax,nzmax,numwfmem) wind components in x-direction [m/s]*
78	! vv(0:nxmax,0:nymax,nzmax,numwfmem) wind components in y-direction [m/s]*
79	! ww(0:nxmax,0:nymax,nzmax,numwfmem) wind components in z-direction *
80	! [deltaeta/s] *
81	! tt(0:nxmax,0:nymax,nzmax,numwfmem) temperature [K] *
82	! pv(0:nxmax,0:nymax,nzmax,numwfmem) potential voriticity (pvu) *
83	! ps(0:nxmax,0:nymax,numwfmem) surface pressure [Pa] *
84	! *
85	!*****************************************************************************
86
87	use par_mod
88	use com_mod
89	use cmapf_mod, only: cc2gll
90
91	implicit none
92
93	real,intent(in),dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: uuh,vvh,pvh
94	real,intent(in),dimension(0:nxmax-1,0:nymax-1,nwzmax) :: wwh
95
96	real,dimension(0:nxmax-1,0:nymax-1,nuvzmax) :: rhoh,uvzlev,wzlev
97	real,dimension(0:nxmax-1,0:nymax-1,nzmax) :: pinmconv
98	!> for reference profile (PS)
99	real :: tvoldref, poldref, pintref, psmean, psstd
100	integer :: ixyref(2)
101	integer, parameter :: ioldref = 1 ! PS 2018-06-13: set to 0 if you
102	!! want old method of searching reference location for the vertical grid
103	!! 1 for new method (options for other methods 2,... in the future)
104
105	real,dimension(0:nymax-1) :: cosf
106	real,dimension(0:nxmax-1,0:nymax-1) :: tvold,pold,pint,tv
107	!! automatic arrays introduced in v10 by ?? to achieve better loop order (PS)
108
109	integer,dimension(0:nxmax-1,0:nymax-1) :: rain_cloud_above,idx
110
111	integer :: ix,jy,kz,iz,n,kmin,ix1,jy1,ixp,jyp,ixref,jyref,kz_inv
112	real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec
113	real :: ew,dz1,dz2,dz
114	real :: xlon,ylat,xlonr,dzdx,dzdy
115	real :: dzdx1,dzdx2,dzdy1,dzdy2
116	real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
117	real,parameter :: const=r_air/ga
118	real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagnostics
119
120	logical, save :: init = .true. ! PS 2018-06-13: add back save attribute
121
122
123	!ZHG SEP 2014 tests
124	! integer :: cloud_ver,cloud_min, cloud_max
125	! integer ::teller(5), convpteller=0, lspteller=0
126	! real :: cloud_col_wat, cloud_water
127	!ZHG 2015 temporary variables for testing
128	! real :: rcw(0:nxmax-1,0:nymax-1)
129	! real :: rpc(0:nxmax-1,0:nymax-1)
130	! character(len=60) :: zhgpath='/xnilu_wrk/flex_wrk/zhg/'
131	! character(len=60) :: fnameA,fnameB,fnameC,fnameD,fnameE,fnameF,fnameG,fnameH
132	! CHARACTER(LEN=3) :: aspec
133	! integer :: virr=0
134	!real :: tot_cloud_h
135	!real :: dbg_height(nzmax)
136	!ZHG
137
138	!*************************************************************************
139	! If verttransform is called the first time, initialize heights of the *
140	! z levels in meter. The heights are the heights of model levels, where *
141	! u,v,T and qv are given, and of the interfaces, where w is given. So, *
142	! the vertical rESOlution in the z system is doubled. As reference point,*
143	! the lower left corner of the grid is used. *
144	! Unlike in the eta system, no difference between heights for u,v and *
145	! heights for w exists. *
146	!*************************************************************************
147
148	!ESO measure CPU time
149	! call mpif_mtime('verttransform',0)
150
151	if (init) then
152
153	!! Search for a point with high surface pressure (i.e. not above significant
154	!! topography) Then, use this point to construct a reference z profile,
155	!! to be used at all times
156	! *****************************************************************************
157
158	if (ioldref .eq. 0) then ! old reference grid point
159	do jy=0,nymin1
160	do ix=0,nxmin1
161	if (ps(ix,jy,1,n).gt.1000.e2) then
162	ixref=ix
163	jyref=jy
164	goto 3
165	endif
166	end do
167	end do
168	3 continue
169
170	! print*,'oldheights at' ,ixref,jyref,ps(ixref,jyref,1,n)
171	else ! new reference grid point
172	! PS: the old version fails if the pressure is <=1000 hPa in the whole
173	! domain. Let us find a good replacement, not just a quick fix.
174	! Search point near to mean pressure after excluding mountains
175
176	psmean = sum( ps(:,:,1,n) ) / (nx*ny)
177	print*,'height: fg psmean',psmean
178	psstd = sqrt(sum( (ps(:,:,1,n) - psmean)*2 ) / (nxny))
179
180	!> new psmean using only points within $\plusminus\sigma$
181	! psmean = sum( ps(:,:,1,n), abs(ps(:,:,1,n)-psmean) < psstd ) / &
182	! count(abs(ps(:,:,1,n)-psmean) < psstd)
183
184	!> new psmean using only points with $p\gt \overbar{p}+\sigma_p$
185	!> (reject mountains, accept valleys)
186	psmean = sum( ps(:,:,1,n), ps(:,:,1,n) > psmean - psstd ) / &
187	count(ps(:,:,1,n) > psmean - psstd)
188	! print*,'height: std, new psmean',psstd,psmean
189	ixyref = minloc( abs( ps(:,:,1,n) - psmean ) )
190	ixref = ixyref(1)
191	jyref = ixyref(2)
192	! print*,'newheights at' ,ixref,jyref,ps(ixref,jyref,1,n)
193	endif
194
195	tvoldref=tt2(ixref,jyref,1,n)* &
196	( 1. + 0.378*ew(td2(ixref,jyref,1,n) ) / ps(ixref,jyref,1,n))
197	poldref=ps(ixref,jyref,1,n)
198	height(1)=0.
199
200	do kz=2,nuvz
201
202	pintref = akz(kz)+bkz(kz)*ps(ixref,jyref,1,n)
203	tv = tth(ixref,jyref,kz,n)(1.+0.608qvh(ixref,jyref,kz,n))
204
205	if (abs(tv(ixref,jyref)-tvoldref) .gt. 0.2) then
206	height(kz) = height(kz-1) + &
207	constlog( poldref/pintref ) &
208	( tv(ixref,jyref) - tvoldref ) / log( tv(ixref,jyref) / tvoldref )
209	else
210	height(kz) = height(kz-1) + &
211	constlog( poldref/pintref ) tv(ixref,jyref)
212	endif
213
214	tvoldref = tv(ixref,jyref)
215	poldref = pintref
216	print*,'height=',kz,height(kz),tvoldref,poldref
217
218	end do
219
220
221	! Determine highest levels that can be within PBL
222	!************************************************
223
224	do kz=1,nz
225	if (height(kz).gt.hmixmax) then
226	nmixz=kz
227	goto 9
228	endif
229	end do
230	9 continue
231
232	! Do not repeat initialization of the Cartesian z grid
233	!*****************************************************
234	init=.false.
235
236	endif ! init block (vertical grid construction)
237
238
239	! Loop over the whole grid
240	!*************************
241
242	tvold(:,:)=tt2(:,:,1,n) * (1.+0.378*ew( td2(:,:,1,n) ) / ps(:,:,1,n))
243	pold=ps(:,:,1,n)
244	uvzlev(:,:,1)=0.
245	wzlev(:,:,1)=0.
246	rhoh(:,:,1)=pold/(r_air*tvold)
247
248
249	! Compute heights of eta levels
250	!******************************
251
252	do kz=2,nuvz
253	pint(:,:)=akz(kz)+bkz(kz)*ps(:,:,1,n)
254	tv(:,:)=tth(:,:,kz,n)(1.+0.608qvh(:,:,kz,n))
255	rhoh(:,:,kz)=pint(:,:)/(r_air*tv(:,:))
256
257	where (abs(tv(:,:)-tvold(:,:)).gt.0.2)
258	uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+constlog(pold(:,:)/pint(:,:)) &
259	(tv(:,:)-tvold(:,:))/log(tv(:,:)/tvold(:,:))
260	elsewhere
261	uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+constlog(pold(:,:)/pint(:,:))tv(:,:)
262	endwhere
263
264	tvold(:,:)=tv(:,:)
265	pold(:,:)=pint(:,:)
266	end do
267
268
269	do kz=2,nwz-1
270	wzlev(:,:,kz)=(uvzlev(:,:,kz+1)+uvzlev(:,:,kz))/2.
271	end do
272	wzlev(:,:,nwz)=wzlev(:,:,nwz-1)+ &
273	uvzlev(:,:,nuvz)-uvzlev(:,:,nuvz-1)
274
275	! pinmconv=(h2-h1)/(p2-p1)
276
277	pinmconv(:,:,1)=(uvzlev(:,:,2))/ &
278	((aknew(2)+bknew(2)*ps(:,:,1,n))- &
279	(aknew(1)+bknew(1)*ps(:,:,1,n)))
280	do kz=2,nz-1
281	pinmconv(:,:,kz)=(uvzlev(:,:,kz+1)-uvzlev(:,:,kz-1))/ &
282	((aknew(kz+1)+bknew(kz+1)*ps(:,:,1,n))- &
283	(aknew(kz-1)+bknew(kz-1)*ps(:,:,1,n)))
284	end do
285	pinmconv(:,:,nz)=(uvzlev(:,:,nz)-uvzlev(:,:,nz-1))/ &
286	((aknew(nz)+bknew(nz)*ps(:,:,1,n))- &
287	(aknew(nz-1)+bknew(nz-1)*ps(:,:,1,n)))
288
289	! Levels where u,v,t and q are given
290	!***********************************
291
292	uu(:,:,1,n)=uuh(:,:,1)
293	vv(:,:,1,n)=vvh(:,:,1)
294	tt(:,:,1,n)=tth(:,:,1,n)
295	qv(:,:,1,n)=qvh(:,:,1,n)
296	!HG adding the cloud water
297	if (readclouds) then
298	clwc(:,:,1,n)=clwch(:,:,1,n)
299	if (.not.sumclouds) ciwc(:,:,1,n)=ciwch(:,:,1,n)
300	end if
301	!HG
302	pv(:,:,1,n)=pvh(:,:,1)
303	rho(:,:,1,n)=rhoh(:,:,1)
304
305	uu(:,:,nz,n)=uuh(:,:,nuvz)
306	vv(:,:,nz,n)=vvh(:,:,nuvz)
307	tt(:,:,nz,n)=tth(:,:,nuvz,n)
308	qv(:,:,nz,n)=qvh(:,:,nuvz,n)
309	!HG adding the cloud water
310	if (readclouds) then
311	clwc(:,:,nz,n)=clwch(:,:,nuvz,n)
312	if (.not. sumclouds) ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n)
313	end if
314	!HG
315	pv(:,:,nz,n)=pvh(:,:,nuvz)
316	rho(:,:,nz,n)=rhoh(:,:,nuvz)
317
318	kmin=2
319	idx=kmin
320	iz_loop: do iz=2,nz-1
321	do jy=0,nymin1
322	do ix=0,nxmin1
323	if (height(iz).gt.uvzlev(ix,jy,nuvz)) then
324
325	uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
326	vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
327	tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
328	qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
329	!HG adding the cloud water
330	if (readclouds) then
331	clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n)
332	if (.not. sumclouds) ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n)
333	end if
334	!HG
335	pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
336	rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
337
338	else
339
340	kz_loop: do kz=idx(ix,jy),nuvz
341	if (idx(ix,jy) .le. kz .and. height(iz).gt.uvzlev(ix,jy,kz-1) &
342	.and. height(iz).le.uvzlev(ix,jy,kz)) then
343	idx(ix,jy)=kz
344	exit kz_loop
345	endif
346	enddo kz_loop
347
348	endif
349	enddo
350	enddo
351	do jy=0,nymin1
352	do ix=0,nxmin1
353	if (height(iz).le.uvzlev(ix,jy,nuvz)) then
354
355	kz=idx(ix,jy)
356	dz1=height(iz)-uvzlev(ix,jy,kz-1)
357	dz2=uvzlev(ix,jy,kz)-height(iz)
358	dz=dz1+dz2
359	uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)dz2+uuh(ix,jy,kz)dz1)/dz
360	vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)dz2+vvh(ix,jy,kz)dz1)/dz
361	tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
362	+tth(ix,jy,kz,n)*dz1)/dz
363	qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
364	+qvh(ix,jy,kz,n)*dz1)/dz
365	!HG adding the cloud water
366	if (readclouds) then
367	clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)dz2+clwch(ix,jy,kz,n)dz1)/dz
368	if (.not. sumclouds) ciwc(ix,jy,iz,n)= &
369	(ciwch(ix,jy,kz-1,n)dz2+ciwch(ix,jy,kz,n)dz1)/dz
370	end if
371	!HG
372	pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)dz2+pvh(ix,jy,kz)dz1)/dz
373	rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)dz2+rhoh(ix,jy,kz)dz1)/dz
374
375	endif
376	enddo
377	enddo
378	enddo iz_loop
379
380
381	! Levels where w is given
382	!*************************
383
384	ww(:,:,1,n)=wwh(:,:,1)*pinmconv(:,:,1)
385	ww(:,:,nz,n)=wwh(:,:,nwz)*pinmconv(:,:,nz)
386	kmin=2
387	idx=kmin
388	iz_loop2: do iz=2,nz
389	do jy=0,nymin1
390	do ix=0,nxmin1
391	kz_loop2: do kz=idx(ix,jy),nwz
392	if (idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1) &
393	.and. height(iz).le.wzlev(ix,jy,kz)) then
394	idx(ix,jy)=kz
395	exit kz_loop2
396	endif
397	enddo kz_loop2
398	enddo
399	enddo
400	do jy=0,nymin1
401	do ix=0,nxmin1
402	kz=idx(ix,jy)
403	dz1=height(iz)-wzlev(ix,jy,kz-1)
404	dz2=wzlev(ix,jy,kz)-height(iz)
405	dz=dz1+dz2
406	ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)pinmconv(ix,jy,kz-1)dz2 &
407	+ wwh(ix,jy,kz) pinmconv(ix,jy,kz)dz1)/dz
408	enddo
409	enddo
410	enddo iz_loop2
411
412	! Compute density gradients at intermediate levels
413	!*************************************************
414
415	drhodz(:,:,1,n)=(rho(:,:,2,n)-rho(:,:,1,n))/(height(2)-height(1))
416	do kz=2,nz-1
417	drhodz(:,:,kz,n)=(rho(:,:,kz+1,n)-rho(:,:,kz-1,n))/ &
418	(height(kz+1)-height(kz-1))
419	end do
420	drhodz(:,:,nz,n)=drhodz(:,:,nz-1,n)
421
422
423	!****************************************************************
424	! Compute slope of eta levels in windward direction and resulting
425	! vertical wind correction
426	!****************************************************************
427
428	do jy=1,ny-2
429	cosf(jy) = 1. / cos( ( real(jy)dy + ylat0 )pi180 )
430	enddo
431
432	kmin=2
433	idx=kmin
434	iz_loop3: do iz=2,nz-1
435	do jy=1,ny-2
436	do ix=1,nx-2
437
438	kz_loop3: do kz=idx(ix,jy),nz
439	if (idx(ix,jy) .le. kz .and. height(iz).gt.uvzlev(ix,jy,kz-1) &
440	.and. height(iz).le.uvzlev(ix,jy,kz)) then
441	idx(ix,jy)=kz
442	exit kz_loop3
443	endif
444	enddo kz_loop3
445	enddo
446	enddo
447
448	do jy=1,ny-2
449	do ix=1,nx-2
450	kz=idx(ix,jy)
451	dz1=height(iz)-uvzlev(ix,jy,kz-1)
452	dz2=uvzlev(ix,jy,kz)-height(iz)
453	dz=dz1+dz2
454	ix1=ix-1
455	jy1=jy-1
456	ixp=ix+1
457	jyp=jy+1
458
459	dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
460	dzdx2=(uvzlev(ixp,jy,kz) -uvzlev(ix1,jy,kz) )/2.
461	dzdx=(dzdx1dz2+dzdx2dz1)/dz
462
463	dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2.
464	dzdy2=(uvzlev(ix,jyp,kz) -uvzlev(ix,jy1,kz) )/2.
465	dzdy=(dzdy1dz2+dzdy2dz1)/dz
466
467	ww(ix,jy,iz,n)=ww(ix,jy,iz,n) +( dzdxuu(ix,jy,iz,n)dxconst*cosf(jy) &
468	+ dzdyvv(ix,jy,iz,n)dyconst)
469
470	enddo
471	enddo
472
473	enddo iz_loop3
474
475	! If north pole is in the domain, calculate wind velocities in polar
476	! stereographic coordinates
477	!*******************************************************************
478
479	if (nglobal) then
480
481	do iz=1,nz
482	do jy=int(switchnorthg)-2,nymin1
483	ylat=ylat0+real(jy)*dy
484	do ix=0,nxmin1
485	xlon=xlon0+real(ix)*dx
486	call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
487	vv(ix,jy,iz,n),uupol(ix,jy,iz,n), vvpol(ix,jy,iz,n))
488	enddo
489	enddo
490	enddo
491
492
493	do iz=1,nz
494
495	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
496	!
497	! AMS nauffer Nov 18 2004 Added check for case vv=0
498
499	xlon=xlon0+real(nx/2-1)*dx
500	xlonr=xlon*pi/180.
501	ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)2+vv(nx/2-1,nymin1,iz,n)2)
502	if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
503	ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr
504	else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
505	ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,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, vvpolaux)
519
520	jy=nymin1
521	do ix=0,nxmin1
522	uupol(ix,jy,iz,n)=uupolaux
523	vvpol(ix,jy,iz,n)=vvpolaux
524	enddo
525	enddo
526
527	! Fix: Set W (vertical wind) at pole to the zonally averaged W of the next
528	! equator-ward parallel
529
530	do iz=1,nz
531	wdummy=0.
532	jy=ny-2
533	do ix=0,nxmin1
534	wdummy=wdummy+ww(ix,jy,iz,n)
535	enddo
536	wdummy=wdummy/real(nx)
537	jy=nymin1
538	do ix=0,nxmin1
539	ww(ix,jy,iz,n)=wdummy
540	enddo
541	enddo
542
543	endif
544
545
546	! If south pole is in the domain, calculate wind velocities in polar
547	! stereographic coordinates
548	!*******************************************************************
549
550	if (sglobal) then
551
552	do iz=1,nz
553	do jy=0,int(switchsouthg)+3
554	ylat=ylat0+real(jy)*dy
555	do ix=0,nxmin1
556	xlon=xlon0+real(ix)*dx
557	call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
558	vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n))
559	enddo
560	enddo
561	enddo
562
563	do iz=1,nz
564
565	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
566	!
567	! AMS nauffer Nov 18 2004 Added check for case vv=0
568
569	xlon=xlon0+real(nx/2-1)*dx
570	xlonr=xlon*pi/180.
571	ffpol=sqrt(uu(nx/2-1,0,iz,n)2+vv(nx/2-1,0,iz,n)2)
572	if (vv(nx/2-1,0,iz,n).lt.0.) then
573	ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
574	else if (vv(nx/2-1,0,iz,n).gt.0.) then
575	ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr
576	else
577	ddpol=pi/2.-xlonr
578	endif
579	if (ddpol.lt.0.) ddpol=2.0*pi+ddpol
580	if (ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
581
582	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
583	xlon=180.0
584	xlonr=xlon*pi/180.
585	ylat=-90.0
586	uuaux=+ffpol*sin(xlonr-ddpol)
587	vvaux=-ffpol*cos(xlonr-ddpol)
588	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux)
589
590	jy=0
591	do ix=0,nxmin1
592	uupol(ix,jy,iz,n)=uupolaux
593	vvpol(ix,jy,iz,n)=vvpolaux
594	enddo
595	enddo
596
597	! Fix: Set W (vertical wind) at pole to the zonally averaged W of the next
598	! equator-ward parallel
599
600	do iz=1,nz
601	wdummy=0.
602	jy=1
603	do ix=0,nxmin1
604	wdummy=wdummy+ww(ix,jy,iz,n)
605	enddo
606	wdummy=wdummy/real(nx)
607	jy=0
608	do ix=0,nxmin1
609	ww(ix,jy,iz,n)=wdummy
610	enddo
611	enddo
612	endif
613
614
615	!******************************************************************************
616	if (readclouds) then ! HG METHOD
617
618	! Loops over all grid cells vertically and construct the 3D matrix for clouds
619	! Cloud top and cloud bottom grid cells are assigned as well as the total column
620	! cloud water. For precipitating columns, the type and whether it is in or below
621	! cloud scavenging are assigned with numbers 2-5 (following the old metod).
622	! A distinction is made between lsp and convp though they are treated the equally
623	! with regard to scavenging. Also, clouds that are not precipitating are defined which
624	! may be used in the future for cloud processing by non-precipitating-clouds.
625	!*******************************************************************************
626
627	!PS write(,) 'Global ECMWF fields: using cloud water'
628	clw(:,:,:,n)=0.0
629	ctwc(:,:,n)=0.0
630	clouds(:,:,:,n)=0
631	! If water/ice are read separately into clwc and ciwc, store sum in clwc
632	if (.not. sumclouds) then
633	clwc(:,:,:,n) = clwc(:,:,:,n) + ciwc(:,:,:,n)
634	endif
635	do jy=0,nymin1
636	do ix=0,nxmin1
637	lsp=lsprec(ix,jy,1,n)
638	convp=convprec(ix,jy,1,n)
639	prec=lsp+convp ! Note PS: prec is not used currently
640	! tot_cloud_h=0
641	! Find clouds in the vertical
642	!! Note PS: bad loop order.
643	do kz=1, nz-1 !go from top to bottom
644	if (clwc(ix,jy,kz,n).gt.0) then
645	! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
646	clw(ix,jy,kz,n)=(clwc(ix,jy,kz,n)rho(ix,jy,kz,n)) &
647	(height(kz+1)-height(kz))
648	! tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
649
650	! icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n) ! Column cloud water [m3/m3]
651	ctwc(ix,jy,n) = ctwc(ix,jy,n)+clw(ix,jy,kz,n)
652	! icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz)) ! Cloud BOT height stats [m]
653	cloudh_min=min(height(kz+1),height(kz))
654	!ZHG 2015 extra for testing
655	! clh(ix,jy,kz,n)=height(kz+1)-height(kz)
656	! icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
657	! icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz)) ! Cloud TOP height [m]
658	!ZHG
659	endif
660	end do
661
662	! If Precipitation. Define removal type in the vertical
663	if (lsp.gt.0.01 .or. convp.gt.0.01) then ! cloud and precipitation
664	!! Note PS: such hardcoded limits would better be parameters defined in par_mod
665
666	do kz=nz,1,-1 !go Bottom up!
667	!! Note PS: bad loop order
668	if (clw(ix,jy,kz,n).gt. 0) then ! is in cloud
669	cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+height(kz)-height(kz-1)
670	clouds(ix,jy,kz,n)=1 ! is a cloud
671	if (lsp.ge.convp) then
672	clouds(ix,jy,kz,n)=3 ! lsp in-cloud
673	else
674	clouds(ix,jy,kz,n)=2 ! convp in-cloud
675	endif ! convective or large scale
676	elseif( clw(ix,jy,kz,n).le.0 .and. cloudh_min.ge.height(kz)) then ! is below cloud
677	if (lsp.ge.convp) then
678	clouds(ix,jy,kz,n)=5 ! lsp dominated washout
679	else
680	clouds(ix,jy,kz,n)=4 ! convp dominated washout
681	endif ! convective or large scale
682	endif
683
684	if (height(kz).ge. 19000) then ! set a max height for removal
685	!! Note PS: such hardcoded limits would better be parameters defined in par_mod
686	clouds(ix,jy,kz,n)=0
687	endif ! clw>0
688	enddo ! kz
689	endif ! precipitation
690	end do
691	end do
692
693	! ESO: copy the relevant data to clw4 to reduce amount of communicated data for MPI
694	! ctwc(:,:,n) = icloud_stats(:,:,4,n)
695
696	!**************************************************************************
697	else ! use old definitions
698	!**************************************************************************
699
700	! create a cloud and rainout/washout field, clouds occur where rh>80%
701	! total cloudheight is stored at level 0
702
703	!PS write(,) 'Global fields: using cloud water from Parameterization'
704	do jy=0,nymin1
705	do ix=0,nxmin1
706	rain_cloud_above(ix,jy)=0
707	lsp=lsprec(ix,jy,1,n)
708	convp=convprec(ix,jy,1,n)
709	cloudsh(ix,jy,n)=0
710	do kz_inv=1,nz-1
711	!! Note PS: bad loop order.
712	kz=nz-kz_inv+1
713	pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
714	rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
715	clouds(ix,jy,kz,n)=0
716
717	if (rh .gt. 0.8) then ! in cloud
718	!! Note PS: such hardcoded limits would better be parameters defined in par_mod
719
720	if (lsp.gt.0.01 .or. convp.gt.0.01) then ! cloud and precipitation
721	!! Note PS: such hardcoded limits would better be parameters defined in par_mod
722	rain_cloud_above(ix,jy)=1
723	cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
724	height(kz)-height(kz-1)
725	if (lsp.ge.convp) then
726	clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
727	else
728	clouds(ix,jy,kz,n)=2 ! convp dominated rainout
729	endif
730	else ! no precipitation
731	clouds(ix,jy,kz,n)=1 ! cloud
732	endif
733
734	else ! no cloud
735
736	if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
737	if (lsp.ge.convp) then
738	clouds(ix,jy,kz,n)=5 ! lsp dominated washout
739	else
740	clouds(ix,jy,kz,n)=4 ! convp dominated washout
741	endif
742	endif
743
744	endif
745	end do
746
747	end do
748	end do
749	endif ! END OLD METHOD
750
751
752	!******* TEST *************
753	! WRITE OUT SOME TEST VARIABLES
754	!******* TEST ********'
755	!teller(:)=0
756	!virr=virr+1
757	!WRITE(aspec, '(i3.3)'), virr
758
759	!if (readclouds) then
760	!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement.txt'
761	!else
762	!fnameH=trim(zhgpath)//trim(aspec)//'Vertical_placement_old.txt'
763	!endif
764	!
765	!OPEN(UNIT=118, FILE=fnameH,FORM='FORMATTED',STATUS = 'UNKNOWN')
766	!do kz_inv=1,nz-1
767	! kz=nz-kz_inv+1
768	! !kz=91
769	! do jy=0,nymin1
770	! do ix=0,nxmin1
771	! if (clouds(ix,jy,kz,n).eq.1) teller(1)=teller(1)+1 ! no precipitation cloud
772	! if (clouds(ix,jy,kz,n).eq.2) teller(2)=teller(2)+1 ! convp dominated rainout
773	! if (clouds(ix,jy,kz,n).eq.3) teller(3)=teller(3)+1 ! lsp dominated rainout
774	! if (clouds(ix,jy,kz,n).eq.4) teller(4)=teller(4)+1 ! convp dominated washout
775	! if (clouds(ix,jy,kz,n).eq.5) teller(5)=teller(5)+1 ! lsp dominated washout
776	!
777	! ! write(,) height(kz),teller
778	! end do
779	! end do
780	! write(118,*) height(kz),teller
781	! teller(:)=0
782	!end do
783	!teller(:)=0
784	!write(,) teller
785	!write(,) aspec
786	!
787	!fnameA=trim(zhgpath)//trim(aspec)//'cloudV.txt'
788	!fnameB=trim(zhgpath)//trim(aspec)//'cloudT.txt'
789	!fnameC=trim(zhgpath)//trim(aspec)//'cloudB.txt'
790	!fnameD=trim(zhgpath)//trim(aspec)//'cloudW.txt'
791	!fnameE=trim(zhgpath)//trim(aspec)//'old_cloudV.txt'
792	!fnameF=trim(zhgpath)//trim(aspec)//'lsp.txt'
793	!fnameG=trim(zhgpath)//trim(aspec)//'convp.txt'
794	!if (readclouds) then
795	!OPEN(UNIT=111, FILE=fnameA,FORM='FORMATTED',STATUS = 'UNKNOWN')
796	!OPEN(UNIT=112, FILE=fnameB,FORM='FORMATTED',STATUS = 'UNKNOWN')
797	!OPEN(UNIT=113, FILE=fnameC,FORM='FORMATTED',STATUS = 'UNKNOWN')
798	!OPEN(UNIT=114, FILE=fnameD,FORM='FORMATTED',STATUS = 'UNKNOWN')
799	!else
800	!OPEN(UNIT=115, FILE=fnameE,FORM='FORMATTED',STATUS = 'UNKNOWN')
801	!OPEN(UNIT=116, FILE=fnameF,FORM='FORMATTED',STATUS = 'UNKNOWN')
802	!OPEN(UNIT=117, FILE=fnameG,FORM='FORMATTED',STATUS = 'UNKNOWN')
803	!endif
804	!
805	!do ix=0,nxmin1
806	!if (readclouds) then
807	!write(111,*) (icloud_stats(ix,jy,1,n),jy=0,nymin1)
808	!write(112,*) (icloud_stats(ix,jy,2,n),jy=0,nymin1)
809	!write(113,*) (icloud_stats(ix,jy,3,n),jy=0,nymin1)
810	!write(114,*) (icloud_stats(ix,jy,4,n),jy=0,nymin1)
811	!else
812	!write(115,*) (cloudsh(ix,jy,n),jy=0,nymin1) !integer
813	!write(116,*) (lsprec(ix,jy,1,n),jy=0,nymin1) !7.83691406E-02
814	!write(117,*) (convprec(ix,jy,1,n),jy=0,nymin1) !5.38330078E-02
815	!endif
816	!end do
817	!
818	!if (readclouds) then
819	!CLOSE(111)
820	!CLOSE(112)
821	!CLOSE(113)
822	!CLOSE(114)
823	!else
824	!CLOSE(115)
825	!CLOSE(116)
826	!CLOSE(117)
827	!endif
828	!
829	!END ******* TEST ************* END
830	! WRITE OUT SOME TEST VARIABLES
831	!END ******* TEST ************* END
832
833
834	! PS 2012
835	! lsp=lsprec(ix,jy,1,n)
836	! convp=convprec(ix,jy,1,n)
837	! prec=lsp+convp
838	! if (lsp.gt.convp) then ! prectype='lsp'
839	! lconvectprec = .false.
840	! else ! prectype='cp'
841	! lconvectprec = .true.
842	! endif
843	! else ! windfields does not contain cloud data
844	! rhmin = 0.90 ! standard condition for presence of clouds
845	!PS note that original by Sabine Eckhart was 80%
846	!PS however, for T<-20 C we consider saturation over ice
847	!PS so I think 90% should be enough
848	! icloudbot(ix,jy,n)=icmv
849	! icloudtop=icmv ! this is just a local variable
850	!98 do kz=1,nz
851	! pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
852	! rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
853	!ps if (prec.gt.0.01) print*,'relhum',prec,kz,rh,height(kz)
854	! if (rh .gt. rhmin) then
855	! if (icloudbot(ix,jy,n) .eq. icmv) then
856	! icloudbot(ix,jy,n)=nint(height(kz))
857	! endif
858	! icloudtop=nint(height(kz)) ! use int to save memory
859	! endif
860	! end do
861	!PS try to get a cloud thicker than 50 m
862	!PS if there is at least .01 mm/h - changed to 0.002 and put into
863	!PS parameter precpmin
864	! if ((icloudbot(ix,jy,n) .eq. icmv .or. &
865	! icloudtop-icloudbot(ix,jy,n) .lt. 50) .and. &
866	! prec .gt. precmin) then
867	! rhmin = rhmin - 0.05
868	! if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum.
869	! end if
870
871	!PS is based on looking at a limited set of comparison data
872	! if (lconvectprec .and. icloudtop .lt. 6000 .and. &
873	! prec .gt. precmin) then
874	!
875	! if (convp .lt. 0.1) then
876	! icloudbot(ix,jy,n) = 500
877	! icloudtop = 8000
878	! else
879	! icloudbot(ix,jy,n) = 0
880	! icloudtop = 10000
881	! endif
882	! endif
883	! if (icloudtop .ne. icmv) then
884	! icloudthck(ix,jy,n) = icloudtop-icloudbot(ix,jy,n)
885	! else
886	! icloudthck(ix,jy,n) = icmv
887	! endif
888	!PS get rid of too thin clouds
889	! if (icloudthck(ix,jy,n) .lt. 50) then
890	! icloudbot(ix,jy,n)=icmv
891	! icloudthck(ix,jy,n)=icmv
892	! endif
893	!HG__________________________________
894	! rcw(ix,jy)=2E-7prec*0.36
895	! rpc(ix,jy)=prec
896	!HG end______________________________
897
898	! endif !HG read clouds
899
900
901
902
903	!ESO measure CPU time
904	! call mpif_mtime('verttransform',1)
905
906	!ESO print out the same measure as in Leo's routine
907	! write(,) 'verttransform: ', &
908	! sum(tt(:,:,:,n)*tt(:,:,:,n)), &
909	! sum(uu(:,:,:,n)uu(:,:,:,n)),sum(vv(:,:,:,n)vv(:,:,:,n)),&
910	! sum(qv(:,:,:,n)qv(:,:,:,n)),sum(pv(:,:,:,n)pv(:,:,:,n)),&
911	! sum(rho(:,:,:,n)rho(:,:,:,n)),sum(drhodz(:,:,:,n)drhodz(:,:,:,n)),&
912	! sum(ww(:,:,:,n)*ww(:,:,:,n)), &
913	! sum(clouds(:,:,:,n)), sum(cloudsh(:,:,n)),sum(idx),sum(pinmconv)
914	end subroutine verttransform_ecmwf
915

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

Download in other formats:

Original Format