Context Navigation

verttransform_ecmwf.f90

10.4.1_peseibugfixes+enhancementsrelease-10.4.1scaling-bug

Last change on this file was fe118c0, checked in by pesei <petra.seibert at univie.ac.at>, 4 years ago
Correct makefile w.r.t. netcdf yes/no and grib_api dependencies Cosmetic correction on verttransform_ecmwf.f90
Property mode set to `100644`
File size: 29.3 KB

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

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

Download in other formats:

Original Format