Context Navigation

verttransform.f90 @ 54cbd6c

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

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

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

Download in other formats:

Original Format