Context Navigation

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

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

Download in other formats:

Original Format