Context Navigation

verttransform_ecmwf.f90 @ 1228ef7

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

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

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

Download in other formats:

Original Format