Context Navigation

verttransform_omega.f90 @ 5f2e8f6

flexpart-noresm

Last change on this file since 5f2e8f6 was 5f2e8f6, checked in by Ignacio Pisso <Ignacio.Pisso@…>, 8 years ago
new flexpart noresm code in src
Property mode set to `100755`
File size: 24.2 KB

Line
1	!**********************************************************************
2	! Copyright 2016 *
3	! Andreas Stohl, Massimo Cassiani, Petra Seibert, A. Frank, *
4	! Gerhard Wotawa, Caroline Forster, Sabine Eckhardt, John Burkhart, *
5	! Harald Sodemann, Ignacio Pisso *
6	! *
7	! This file is part of FLEXPART-NorESM *
8	! *
9	! FLEXPART-NorESM is free software: you can redistribute it *
10	! and/or modify *
11	! it under the terms of the GNU General Public License as published by*
12	! the Free Software Foundation, either version 3 of the License, or *
13	! (at your option) any later version. *
14	! *
15	! FLEXPART-NorESM is distributed in the hope that it will be useful, *
16	! but WITHOUT ANY WARRANTY; without even the implied warranty of *
17	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
18	! GNU General Public License for more details. *
19	! *
20	! You should have received a copy of the GNU General Public License *
21	! along with FLEXPART-NorESM. *
22	! If not, see <http://www.gnu.org/licenses/>. *
23	!**********************************************************************
24
25	subroutine verttransform_omega(n,uuh,vvh,wwh,pvh,itime,indj)
26	! i i i i i i i
27	!*****************************************************************************
28	! *
29	! This subroutine transforms temperature, dew point temperature and *
30	! wind components from eta to meter coordinates. *
31	! The vertical wind component is transformed from Pa/s to m/s using *
32	! the conversion factor pinmconv. *
33	! In addition, this routine calculates vertical density gradients *
34	! needed for the parameterization of the turbulent velocities. *
35	! *
36	! Author: A. Stohl, G. Wotawa *
37	! *
38	! 12 August 1996 *
39	! Update: 16 January 1998 *
40	! *
41	! Major update: 17 February 1999 *
42	! by G. Wotawa *
43	! *
44	! - Vertical levels for u, v and w are put together *
45	! - Slope correction for vertical velocity: Modification of calculation *
46	! procedure *
47	!*****************************************************************************
48	! Modified by M. Cassiani 2106 to trasnform omega in NorESM in w and *
49	! interpolate on FLEXPART levels, afterwards put into *
50	! terrain following coordinate *
51	! *
52	!*****************************************************************************
53	! Changes, Bernd C. Krueger, Feb. 2001:
54	! Variables tth and qvh (on eta coordinates) from common block
55	!*****************************************************************************
56	! Sabine Eckhardt, March 2007
57	! added the variable cloud for use with scavenging - descr. in com_mod
58	!*****************************************************************************
59	! *
60	! Variables: *
61	! nx,ny,nz field dimensions in x,y and z direction *
62	! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition *
63	! uu(0:nxmax,0:nymax,nzmax,2) wind components in x-direction [m/s] *
64	! vv(0:nxmax,0:nymax,nzmax,2) wind components in y-direction [m/s] *
65	! ww(0:nxmax,0:nymax,nzmax,2) wind components in z-direction [deltaeta/s]*
66	! tt(0:nxmax,0:nymax,nzmax,2) temperature [K] *
67	! pv(0:nxmax,0:nymax,nzmax,2) potential voriticity (pvu) *
68	! ps(0:nxmax,0:nymax,2) 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	integer, save :: indice(3)=0,helpint=0,counter=0,izp,iz1,tempo(0:3)=0
79	integer :: indj
80	integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
81	integer :: rain_cloud_above,kz_inv
82	real :: f_qvsat,pressure
83	real :: rh,lsp,convp
84	real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
85	real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
86	real :: xlon,ylat,xlonr,dzdx,dzdy
87	real :: dzdx1,dzdx2,dzdy1,dzdy2
88	real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
89	real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
90	real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
91	real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
92	real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
93	real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
94
95
96	real,parameter :: const=r_air/ga
97
98	logical :: init = .true.
99	!below added by mc for test reason
100	real latitudine, longitudine, quotap, quotaz
101	character :: adate8,atime6
102	real(kind=dp) :: jul
103	integer :: itime,i,j,jjjjmmdd,ihmmss
104	integer :: skyp_check_w=1 !,timeintervalins=86400 !!this flag test w, timeintervalins must be an integer multiple of the input tiem interval e.g. 3h
105	character*7 :: stringwftime ! added for testing by mc
106
107	!C------------------ TEST FILE TO BE DELETED IN FINAL VERSION ------------
108
109	if (skyp_check_w.eq.0) then
110	write(stringwftime,'(I7.7)')wftime(indj)
111	open(79,file='..\windtrans\w_from_omega__'//stringwftime//'.dat') !
112	! write(stringwftime,'(I7.7)')wftime(indj)
113	!open(79,file='..\windtrans_omega\w_fromomega'//stringwftime//'.dat') !
114	end if
115	! c-------------------------------------------------------------------------
116	! Determine current calendar date
117	!**********************************************************
118	jul=bdate+real(itime,kind=dp)/86400._dp
119	call caldate(jul,jjjjmmdd,ihmmss)
120	write(adate,'(i8.8)') jjjjmmdd
121	write(atime,'(i6.6)') ihmmss
122	! Open output file and write the output
123	!**************************************
124	! open(137,file=path(2)(1:length(2))//'testveloc_'//adate// & ! for testing u v interpolation against original data : by mc
125	! atime,form='formatted')
126	!
127	!100 format(7e20.10)
128	!******* above added for test reSAON BY MC ***************
129
130	!*************************************************************************
131	! If verttransform is called the first time, initialize heights of the *
132	! z levels in meter. The heights are the heights of model levels, where *
133	! u,v,T and qv are given, and of the interfaces, where w is given. So, *
134	! the vertical resolution in the z system is doubled. As reference point,*
135	! the lower left corner of the grid is used. *
136	! Unlike in the eta system, no difference between heights for u,v and *
137	! heights for w exists. *
138	!*************************************************************************
139
140
141	!do kz=1,nuvz
142	! write (,) 'akz: ',akz(kz),'bkz',bkz(kz)
143	!end do
144
145	if (init) then
146
147	! Search for a point with high surface pressure (i.e. not above significant topography)
148	! Then, use this point to construct a reference z profile, to be used at all times
149	!*****************************************************************************
150
151	do jy=0,nymin1
152	do ix=0,nxmin1
153	if (ps(ix,jy,1,n).gt.100000.) then
154	ixm=ix
155	jym=jy
156	goto 3
157	endif
158	end do
159	end do
160	3 continue
161
162
163	tvold=tt2(ixm,jym,1,n)(1.+0.378ew(td2(ixm,jym,1,n))/ &
164	ps(ixm,jym,1,n))
165	pold=ps(ixm,jym,1,n)
166	height(1)=0.
167
168	do kz=2,nuvz
169	pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
170	tv=tth(ixm,jym,kz,n)(1.+0.608qvh(ixm,jym,kz,n))
171
172
173	! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
174	! upon the transformation to z levels. In order to save computer memory, this is
175	! not done anymore in the standard version. However, this option can still be
176	! switched on by replacing the following lines with those below, that are
177	! currently commented out.
178	! Note that two more changes are necessary in this subroutine below.
179	! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
180	!*****************************************************************************
181
182	if (abs(tv-tvold).gt.0.2) then
183	height(kz)= &
184	height(kz-1)+constlog(pold/pint) &
185	(tv-tvold)/log(tv/tvold)
186	else
187	height(kz)=height(kz-1)+ &
188	constlog(pold/pint)tv
189	endif
190
191	! Switch on following lines to use doubled vertical resolution
192	!*************************************************************
193	! if (abs(tv-tvold).gt.0.2) then
194	! height((kz-1)*2)=
195	! + height(max((kz-2)2,1))+constlog(pold/pint)*
196	! + (tv-tvold)/log(tv/tvold)
197	! else
198	! height((kz-1)2)=height(max((kz-2)2,1))+
199	! + constlog(pold/pint)tv
200	! endif
201	! End doubled vertical resolution
202
203	tvold=tv
204	pold=pint
205	end do
206
207
208	! Switch on following lines to use doubled vertical resolution
209	!*************************************************************
210	! do 7 kz=3,nz-1,2
211	! height(kz)=0.5*(height(kz-1)+height(kz+1))
212	! height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
213	! End doubled vertical resolution
214
215
216	! Determine highest levels that can be within PBL
217	!************************************************
218
219	do kz=1,nz
220	if (height(kz).gt.hmixmax) then
221	nmixz=kz
222	goto 9
223	endif
224	end do
225	9 continue
226
227	! Do not repeat initialization of the Cartesian z grid
228	!*****************************************************
229
230	init=.false.
231
232	endif !on init
233
234
235	! Loop over the whole grid
236	!*************************
237
238	do jy=0,nymin1
239	do ix=0,nxmin1
240	latitudine= xlon0+(ix)*dx !added for test reason by mc
241	longitudine=ylat0+(jy)*dy !added for test reason by mc
242
243	tvold=tt2(ix,jy,1,n)(1.+0.378ew(td2(ix,jy,1,n))/ &
244	ps(ix,jy,1,n))
245	pold=ps(ix,jy,1,n)
246	uvzlev(1)=0.
247	wzlev(1)=0.
248	rhoh(1)=pold/(r_air*tvold)
249	!rhoperu(ix,jy,1)=uuh(ix,jy,1)*rhoh(1) !aaded for testing by mc 15-11-2013
250	!rhoperv(ix,jy,1)=vvh(ix,jy,1)*rhoh(1) !aaded for testing by mc 15-11-2013
251
252	! Compute heights of eta levels
253	!******************************
254
255	do kz=2,nuvz
256	pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
257	tv=tth(ix,jy,kz,n)(1.+0.608qvh(ix,jy,kz,n))
258	rhoh(kz)=pint/(r_air*tv)
259	!rhoperu(ix,jy,kz)=uuh(ix,jy,kz)*rhoh(kz) !aaded for testing by mc 15-11-2013
260	!rhoperv(ix,jy,kz)=vvh(ix,jy,kz)*rhoh(kz) !aaded for testing by mc 15-11-2013
261	if (abs(tv-tvold).gt.0.2) then
262	uvzlev(kz)=uvzlev(kz-1)+constlog(pold/pint) &
263	(tv-tvold)/log(tv/tvold)
264	else
265	uvzlev(kz)=uvzlev(kz-1)+constlog(pold/pint)tv
266	endif
267
268	tvold=tv
269	pold=pint
270	end do
271
272
273	do kz=2,netadot-1 !--NOTE: netadot is used in NORESM in ECMWF it was nwz. netadot is used since omega in Pa/s
274	!--is co-located with u and etadot levels (akm, bkm) are not co-located and therefore used for computing pressur gradients
275	wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2. !high of etadot levels
276	end do
277	wzlev(netadot)=wzlev(netadot-1)+ &
278	uvzlev(nuvz)-uvzlev(nuvz-1)
279
280	uvwzlev(ix,jy,1)=0.0
281	do kz=2,nuvz
282	uvwzlev(ix,jy,kz)=uvzlev(kz)
283	end do
284
285	! doubling of vertical resolution is not implemeted yet in NORESM version --
286	! Switch on following lines to use doubled vertical resolution
287	! Switch off the three lines above.
288	!*************************************************************
289	!22 uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
290	! do 23 kz=2,nwz
291	!23 uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
292	! End doubled vertical resolution
293
294
295	pinmconv(1)=(uvzlev(2)-wzlev(1))/ & !noRESM. note: akm are etadot levels staggered from u,v,omega levels see e.g. CAM3 user guide page 52
296	((akz(2)+bkz(2)*ps(ix,jy,1,n))- & !this is not actually used
297	(akm(1)+bkm(1)*ps(ix,jy,1,n)))
298
299	do kz=2,nz-1
300	pinmconv(kz)=(wzlev(kz)-wzlev(kz-1))/ & !NORESM
301	((akm(kz)+bkm(kz)*ps(ix,jy,1,n))- &
302	(akm(kz-1)+bkm(kz-1)*ps(ix,jy,1,n)))
303	end do
304	pinmconv(nz)=(wzlev(nz)-wzlev(nz-1))/ & !NORESM
305	((akm(nz)+bkm(nz)*ps(ix,jy,1,n))- &
306	(akm(nz-1)+bkm(nz-1)*ps(ix,jy,1,n)))
307
308	! Levels, where u,v,t and q are given
309	!************************************
310	ww(ix,jy,1,n)=wwh(ix,jy,1)/(-rhoh(1)*ga) ! NORESM note that here teh value of zero is assigned to vertical velocity at the ground see below
311	uu(ix,jy,1,n)=uuh(ix,jy,1)
312	vv(ix,jy,1,n)=vvh(ix,jy,1)
313	tt(ix,jy,1,n)=tth(ix,jy,1,n)
314	qv(ix,jy,1,n)=qvh(ix,jy,1,n)
315	pv(ix,jy,1,n)=pvh(ix,jy,1)
316	rho(ix,jy,1,n)=rhoh(1)
317	uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
318	vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
319	ww(ix,jy,nz,n)=wwh(ix,jy,nuvz)/(-rhoh(nuvz)*ga) !NORESM noet that wwh(..,nz) has not been assigned zero here in readwind (zero assginment done in ECMWF)
320	tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
321	qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
322	pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
323	rho(ix,jy,nz,n)=rhoh(nuvz)
324	kmin=2
325	do iz=2,nz-1
326	do kz=kmin,nuvz
327	if(height(iz).gt.uvzlev(nuvz)) then
328	uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
329	vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
330	tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
331	qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
332	pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
333	rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
334	goto 30
335	endif
336	if ((height(iz).gt.uvzlev(kz-1)).and. &
337	(height(iz).le.uvzlev(kz))) then
338	dz1=height(iz)-uvzlev(kz-1)
339	dz2=uvzlev(kz)-height(iz)
340	dz=dz1+dz2
341	uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)dz2+uuh(ix,jy,kz)dz1)/dz
342	vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)dz2+vvh(ix,jy,kz)dz1)/dz
343	if (height(iz).gt.uvzlev(2)) then !NORESM using cartesian OMEGA=-rhogW by mc
344	ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)/(-rhoh(kz-1)ga)dz2 & !NORESM using cartesian OMEGA=-rhogW
345	+wwh(ix,jy,kz)/(-rhoh(kz)ga)dz1)/dz !NORESM
346	else if (height(iz).le.uvzlev(2)) then !NORESM
347	ww(ix,jy,iz,n)=wwh(ix,jy,2)/(-rhoh(2)*ga) !NORESM
348	!else !NORESM
349	!ww(ix,jy,iz,n)=0. !NORESM ww will be assigned below
350	end if
351
352	! write (137,100)1.*n, latitudine, longitudine, height(iz), uu(ix,jy,iz,n), vv(ix,jy,iz,n), ww(ix,jy,iz,n) !ADDDEd FOR TEST REASON BY MC !TO BE REMOMEVED
353
354	tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
355	+tth(ix,jy,kz,n)*dz1)/dz
356	qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
357	+qvh(ix,jy,kz,n)*dz1)/dz
358	pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)dz2+pvh(ix,jy,kz)dz1)/dz
359	rho(ix,jy,iz,n)=(rhoh(kz-1)dz2+rhoh(kz)dz1)/dz
360	kmin=kz
361	goto 30
362	endif
363	end do
364	30 continue
365	end do
366
367
368
369
370	! Compute density gradients at intermediate levels
371	!*************************************************
372
373	drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
374	(height(2)-height(1))
375	do kz=2,nz-1
376	drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
377	(height(kz+1)-height(kz-1))
378	end do
379	drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
380
381	end do !
382	end do ! on jy=0,nymin1
383
384
385
386	!****************************************************************
387	! Correct w velocity for topogrphy
388	!****************************************************************
389
390	do jy=1,ny-2
391	do ix=1,nx-2
392	ix1=ix-1
393	jy1=jy-1
394	ixp=ix+1
395	jyp=jy+1
396	do iz=2,nz-1
397	ui=uu(ix,jy,iz,n)dxconst/cos((real(jy)dy+ylat0)*pi180)
398	vi=vv(ix,jy,iz,n)*dyconst
399	dzdx=(oro(ixp,jy) - oro(ix1,jy))/2.
400	dzdy=(oro(ix,jyp) - oro(ix,jy1))/2.
401	ww(ix,jy,iz,n)=ww(ix,jy,iz,n)-(dzdxui+dzdyvi)
402	end do
403	end do
404	end do
405	!*******************************************************************
406
407	!-------additional part for NORESM use value of w in terrain following coordinate
408	!-------assign w=0 in terrain following coordinate at terrain level and, to fill the value not already set,
409	!------- i.e value for height less then uvwzle(2), interpolate between zero and the first assigned value
410
411	do jy=1,ny-2
412	do ix=1,nx-2
413	do iz=nz,2,-1
414	if (height(iz).lt.uvwzlev(ix,jy,2)) then
415	dz=uvwzlev(ix,jy,2)-height(1)
416	dz1=height(iz)-height(1)
417	dz2=dz-dz1
418	!print *,ww(ix,jy,iz,n)
419	ww(ix,jy,iz,n)=ww(ix,jy,1,n)dz2/dz+ww(ix,jy,iz,n)dz1/dz
420	!print *,ww(ix,jy,iz,n)
421	end if
422	end do
423	end do
424	end do
425
426
427	!C-------------- write some diagnostic by mc ---------------
428	if (skyp_check_w.eq.0) then
429	!if (mod(wftime(indj),timeintervalins).eq.0) then
430	do jy=0,nymin1
431	do ix=0, nxfield-1
432	do kz=1,26
433	write (79,'(3f12.4,4E15.6)')xlon0+dxix,ylat0+dyjy,height(kz),ww(ix,jy,kz,n)
434	end do
435	end do
436	end do
437	!end if
438	end if
439
440
441	! If north pole is in the domain, calculate wind velocities in polar
442	! stereographic coordinates
443	!*******************************************************************
444
445	if (nglobal) then
446	do jy=int(switchnorthg)-2,nymin1
447	ylat=ylat0+real(jy)*dy
448	do ix=0,nxmin1
449	xlon=xlon0+real(ix)*dx
450	do iz=1,nz
451	call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
452	vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
453	vvpol(ix,jy,iz,n))
454	end do
455	end do
456	end do
457
458
459	do iz=1,nz
460
461	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
462	!
463	! AMSnauffer Nov 18 2004 Added check for case vv=0
464	!
465	xlon=xlon0+real(nx/2-1)*dx
466	xlonr=xlon*pi/180.
467	ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
468	vv(nx/2-1,nymin1,iz,n)**2)
469	if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
470	ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
471	vv(nx/2-1,nymin1,iz,n))-xlonr
472	else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
473	ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
474	vv(nx/2-1,nymin1,iz,n))-xlonr
475	else
476	ddpol=pi/2-xlonr
477	endif
478	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
479	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
480
481	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
482	xlon=180.0
483	xlonr=xlon*pi/180.
484	ylat=90.0
485	uuaux=-ffpol*sin(xlonr+ddpol)
486	vvaux=-ffpol*cos(xlonr+ddpol)
487	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
488	vvpolaux)
489	jy=nymin1
490	do ix=0,nxmin1
491	uupol(ix,jy,iz,n)=uupolaux
492	vvpol(ix,jy,iz,n)=vvpolaux
493	end do
494	end do
495
496
497	! Fix: Set W at pole to the zonally averaged W of the next equator-
498	! ward parallel of latitude
499
500	do iz=1,nz
501	wdummy=0.
502	jy=ny-2
503	do ix=0,nxmin1
504	wdummy=wdummy+ww(ix,jy,iz,n)
505	end do
506	wdummy=wdummy/real(nx)
507	jy=nymin1
508	do ix=0,nxmin1
509	ww(ix,jy,iz,n)=wdummy
510	end do
511	end do
512
513	endif
514
515
516	! If south pole is in the domain, calculate wind velocities in polar
517	! stereographic coordinates
518	!*******************************************************************
519
520	if (sglobal) then
521	do jy=0,int(switchsouthg)+3
522	ylat=ylat0+real(jy)*dy
523	do ix=0,nxmin1
524	xlon=xlon0+real(ix)*dx
525	do iz=1,nz
526	call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
527	vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
528	vvpol(ix,jy,iz,n))
529	end do
530	end do
531	end do
532
533	do iz=1,nz
534
535	! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
536	!
537	! AMSnauffer Nov 18 2004 Added check for case vv=0
538	!
539	xlon=xlon0+real(nx/2-1)*dx
540	xlonr=xlon*pi/180.
541	ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
542	vv(nx/2-1,0,iz,n)**2)
543	if (vv(nx/2-1,0,iz,n).lt.0.) then
544	ddpol=atan(uu(nx/2-1,0,iz,n)/ &
545	vv(nx/2-1,0,iz,n))+xlonr
546	else if (vv(nx/2-1,0,iz,n).gt.0.) then
547	ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
548	vv(nx/2-1,0,iz,n))+xlonr
549	else
550	ddpol=pi/2-xlonr
551	endif
552	if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
553	if(ddpol.gt.2.0pi) ddpol=ddpol-2.0pi
554
555	! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
556	xlon=180.0
557	xlonr=xlon*pi/180.
558	ylat=-90.0
559	uuaux=+ffpol*sin(xlonr-ddpol)
560	vvaux=-ffpol*cos(xlonr-ddpol)
561	call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
562	vvpolaux)
563
564	jy=0
565	do ix=0,nxmin1
566	uupol(ix,jy,iz,n)=uupolaux
567	vvpol(ix,jy,iz,n)=vvpolaux
568	end do
569	end do
570
571
572	! Fix: Set W at pole to the zonally averaged W of the next equator-
573	! ward parallel of latitude
574
575	do iz=1,nz
576	wdummy=0.
577	jy=1
578	do ix=0,nxmin1
579	wdummy=wdummy+ww(ix,jy,iz,n)
580	end do
581	wdummy=wdummy/real(nx)
582	jy=0
583	do ix=0,nxmin1
584	ww(ix,jy,iz,n)=wdummy
585	end do
586	end do
587	endif
588
589
590	!write (,) 'initializing clouds, n:',n,nymin1,nxmin1,nz
591	! create a cloud and rainout/washout field, clouds occur where rh>80%
592	! total cloudheight is stored at level 0
593	do jy=0,nymin1
594	do ix=0,nxmin1
595	rain_cloud_above=0
596	lsp=lsprec(ix,jy,1,n)
597	convp=convprec(ix,jy,1,n)
598	cloudsh(ix,jy,n)=0
599	do kz_inv=1,nz-1
600	kz=nz-kz_inv+1
601	pressure=rho(ix,jy,kz,n)r_airtt(ix,jy,kz,n)
602	rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
603	clouds(ix,jy,kz,n)=0
604	if (rh.gt.0.8) then ! in cloud
605	if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
606	rain_cloud_above=1
607	cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
608	height(kz)-height(kz-1)
609	if (lsp.ge.convp) then
610	clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
611	else
612	clouds(ix,jy,kz,n)=2 ! convp dominated rainout
613	endif
614	else ! no precipitation
615	clouds(ix,jy,kz,n)=1 ! cloud
616	endif
617	else ! no cloud
618	if (rain_cloud_above.eq.1) then ! scavenging
619	if (lsp.ge.convp) then
620	clouds(ix,jy,kz,n)=5 ! lsp dominated washout
621	else
622	clouds(ix,jy,kz,n)=4 ! convp dominated washout
623	endif
624	endif
625	endif
626	end do
627	end do
628	end do
629
630	!do 102 kz=1,nuvz
631	!write(an,'(i02)') kz+10
632	!write(,) nuvz,nymin1,nxmin1,'--',an,'--'
633	!open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
634	!do 101 jy=0,nymin1
635	! write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
636	!101 continue
637	! close(4)
638	!102 continue
639
640	! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
641	! do 103 jy=0,nymin1
642	! write (4,*)
643	!+ (height(kz),kz=1,nuvz)
644	!103 continue
645	! close(4)
646
647	!open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
648	! do 104 jy=0,nymin1
649	! write (4,*)
650	!+ (r_airtt(ix,jy,1,n)rho(ix,jy,1,n),ix=1,nxmin1)
651	!104 continue
652	! close(4)
653	end subroutine verttransform_omega

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

Download in other formats:

Original Format