Context Navigation

verttransform.f90 @ b0434e1

10.4.1_peseiGFS_025bugfixes+enhancementsdevrelease-10release-10.4.1scaling-bugunivie

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

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

Download in other formats:

Original Format