Context Navigation

verttransform_ecmwf.f90 @ 3481cc1

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

Last change on this file since 3481cc1 was 3481cc1, checked in by Ignacio Pisso <ip@…>, 4 years ago
move license from headers to a different file
Property mode set to `100644`
File size: 29.1 KB

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

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

Download in other formats:

Original Format