Context Navigation

verttransform_ecmwf.f90 @ adead08

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

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

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

Download in other formats:

Original Format