Context Navigation

verttransform_nests.f90 @ 1a8fbee

univie

Last change on this file since 1a8fbee was 1a8fbee, checked in by pesei <petra seibert at univie ac at>, 6 years ago
vertransform: fix #140 (ref position for grid), verttransform_nest: fix bug with dimensions, both: smaller changes, com_mod: correct comment to nmixz
Property mode set to `100644`
File size: 24.5 KB

Rev	Line
[e200b7a]	1	!**********************************************************************
	2	! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010 *
	3	! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa, *
	4	! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann *
	5	! *
	6	! This file is part of FLEXPART. *
	7	! *
	8	! FLEXPART is free software: you can redistribute it and/or modify *
	9	! it under the terms of the GNU General Public License as published by*
	10	! the Free Software Foundation, either version 3 of the License, or *
	11	! (at your option) any later version. *
	12	! *
	13	! FLEXPART is distributed in the hope that it will be useful, *
	14	! but WITHOUT ANY WARRANTY; without even the implied warranty of *
	15	! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
	16	! GNU General Public License for more details. *
	17	! *
	18	! You should have received a copy of the GNU General Public License *
	19	! along with FLEXPART. If not, see <http://www.gnu.org/licenses/>. *
	20	!**********************************************************************
	21
	22	subroutine verttransform_nests(n,uuhn,vvhn,wwhn,pvhn)
[db712a8]	23	! i i i i i
	24	!*****************************************************************************
	25	! *
	26	! This subroutine transforms temperature, dew point temperature and *
	27	! wind components from eta to meter coordinates. *
	28	! The vertical wind component is transformed from Pa/s to m/s using *
	29	! the conversion factor pinmconv. *
	30	! In addition, this routine calculates vertical density gradients *
	31	! needed for the parameterization of the turbulent velocities. *
	32	! It is similar to verttransform, but makes the transformations for *
	33	! the nested grids. *
	34	! *
	35	! Author: A. Stohl, G. Wotawa *
	36	! *
	37	! 12 August 1996 *
	38	! Update: 16 January 1998 *
	39	! *
[1a8fbee]	40	! *
	41	!*****************************************************************************
	42	! CHANGES *
[db712a8]	43	! Major update: 17 February 1999 *
	44	! by G. Wotawa *
	45	! *
	46	! - Vertical levels for u, v and w are put together *
	47	! - Slope correction for vertical velocity: Modification of calculation *
	48	! procedure *
	49	! *
	50	!*****************************************************************************
	51	! Changes, Bernd C. Krueger, Feb. 2001: (marked "C-cv")
	52	! Variables tthn and qvhn (on eta coordinates) from common block
[1a8fbee]	53	! Sabine Eckhardt, March 2007:
	54	! added the variable cloud for use with scavenging - descr. in com_mod
	55	! *
	56	! Who? When? *
	57	! Unified ECMWF and GFS builds
	58	! Marian Harustak, 12.5.2017
	59	! - Renamed from verttransform to verttransform_ecmwf
	60	! *
[db712a8]	61	! ESO, 2016
	62	! -note that divide-by-zero occurs when nxmaxn,nymaxn etc. are larger than
[1a8fbee]	63	! the actual field dimensions /fixed, PS 2018/
	64	! *
	65	! Don Morton, 2017-05-30:
	66	! modification of a bug in ew. Don Morton (CTBTO project) *
	67	! *
	68	! undocumented modifications by NILU for v10 *
	69	! *
	70	! Petra Seibert, 2018-06-13: *
	71	! - insert proper boundaries for implied loops in array expressions *s
	72	! - minor changes, most of them just cosmetics *
	73	! for details see changelog.txt *
	74	! *
	75	! ****************************************************************************
[db712a8]	76	! *
	77	! Variables: *
	78	! nxn,nyn,nuvz,nwz field dimensions in x,y and z direction *
	79	! uun wind components in x-direction [m/s] *
	80	! vvn wind components in y-direction [m/s] *
	81	! wwn wind components in z-direction [deltaeta/s]*
	82	! ttn temperature [K] *
	83	! pvn potential vorticity (pvu) *
	84	! psn surface pressure [Pa] *
	85	! *
	86	!*****************************************************************************
[e200b7a]	87
	88	use par_mod
	89	use com_mod
	90
	91	implicit none
	92
[1a8fbee]	93	real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax,maxnests) :: &
	94	uuhn,vvhn,pvhn
[db712a8]	95	real,intent(in),dimension(0:nxmaxn-1,0:nymaxn-1,nwzmax,maxnests) :: wwhn
	96
	97	real,dimension(0:nxmaxn-1,0:nymaxn-1,nuvzmax) :: rhohn,uvzlev,wzlev
	98	real,dimension(0:nxmaxn-1,0:nymaxn-1,nzmax) :: pinmconv
[1a8fbee]	99
[db712a8]	100	real,dimension(0:nymaxn-1) :: cosf
[1a8fbee]	101	real,dimension(0:nxmaxn-1,0:nymaxn-1) :: tvold,pold,pint,tv
	102	!! automatic arrays introduced in v10 by ?? to achieve better loop order (PS)
[db712a8]	103
	104	integer,dimension(0:nxmaxn-1,0:nymaxn-1) :: rain_cloud_above, idx
	105
	106	integer :: ix,jy,kz,iz,n,l,kmin,kl,klp,ix1,jy1,ixp,jyp,kz_inv
[1a8fbee]	107	integer :: nxx, nyy ! max of nest where we are working
[db712a8]	108	real :: f_qvsat,pressure,rh,lsp,convp,cloudh_min,prec
	109
	110	real :: ew,dz1,dz2,dz
[e200b7a]	111	real :: dzdx,dzdy
	112	real :: dzdx1,dzdx2,dzdy1,dzdy2
	113	real,parameter :: const=r_air/ga
[db712a8]	114	real :: tot_cloud_h
[e200b7a]	115
[1a8fbee]	116	! real,parameter :: precmin = 0.002 ! minimum prec in mm/h for cloud diagn.
[e200b7a]	117
[db712a8]	118	! Loop over all nests
	119	!********************
[e200b7a]	120
[1a8fbee]	121	l_loop: do l=1,numbnests
	122	nxx=nxn(l)-1 ! shorthand
	123	nyy=nyn(l)-1 ! shorthand
[e200b7a]	124
[1a8fbee]	125	! Loop over the whole grid (partly implicitly
[db712a8]	126	!*************************
[1a8fbee]	127	! initialise the automatic arrays
[e200b7a]	128
[1a8fbee]	129	tvold(0:nxx,0:nyy)=tt2n(0:nxx,0:nyy,1,n,l) * &
	130	(1.+0.378*ew( td2n(0:nxx,0:nyy,1,n,l) ) / psn(0:nxx,0:nyy,1,n,l))
	131	pold(0:nxx,0:nyy)=psn(0:nxx,0:nyy,1,n,l)
	132	uvzlev(0:nxx,0:nyy,1)=0.
	133	wzlev(0:nxx,0:nyy,1)=0.
	134	rhohn(0:nxx,0:nyy,1)=pold(0:nxx,0:nyy)/(r_air*tvold(0:nxx,0:nyy))
[e200b7a]	135
	136
[db712a8]	137	! Compute heights of eta levels
	138	!******************************
[e200b7a]	139
[db712a8]	140	do kz=2,nuvz
[1a8fbee]	141	pint(:,:)=akz(kz)+bkz(kz)*psn(0:nxx,0:nyy,1,n,l)
	142	tv(0:nxx,0:nyy)=tthn(0:nxx,0:nyy,kz,n,l)* &
	143	(1.+0.608*qvhn(0:nxx,0:nyy,kz,n,l))
	144	rhohn(0:nxx,0:nyy,kz)=pint(0:nxx,0:nyy)/(r_air*tv(0:nxx,0:nyy))
	145
	146	where (abs(tv(0:nxx,0:nyy)-tvold(0:nxx,0:nyy)).gt.0.2)
	147	uvzlev(0:nxx,0:nyy,kz)=uvzlev(0:nxx,0:nyy,kz-1) + &
	148	constlog( pold(0:nxx,0:nyy)/pint(0:nxx,0:nyy) ) &
	149	( tv(0:nxx,0:nyy) - tvold(0:nxx,0:nyy) ) / &
	150	log( tv(0:nxx,0:nyy)/tvold(0:nxx,0:nyy) )
[db712a8]	151	elsewhere
[1a8fbee]	152	uvzlev(0:nxx,0:nyy,kz)=uvzlev(0:nxx,0:nyy,kz-1) + &
	153	constlog( pold(0:nxx,0:nyy)/pint(0:nxx,0:nyy) ) &
	154	tv(0:nxx,0:nyy)
[db712a8]	155	endwhere
[e200b7a]	156
[1a8fbee]	157	tvold(0:nxx,0:nyy)=tv(0:nxx,0:nyy)
	158	pold(0:nxx,0:nyy)=pint(0:nxx,0:nyy)
[e200b7a]	159
[db712a8]	160	end do
[e200b7a]	161
[db712a8]	162	do kz=2,nwz-1
[1a8fbee]	163	wzlev(0:nxx,0:nyy,kz)=(uvzlev(0:nxx,0:nyy,kz+1)+uvzlev(0:nxx,0:nyy,kz))/2.
[db712a8]	164	end do
[1a8fbee]	165	wzlev(0:nxx,0:nyy,nwz)=wzlev(0:nxx,0:nyy,nwz-1)+ &
	166	uvzlev(0:nxx,0:nyy,nuvz)-uvzlev(0:nxx,0:nyy,nuvz-1)
[e200b7a]	167
	168
[1a8fbee]	169	pinmconv(0:nxx,0:nyy,1)=(uvzlev(0:nxx,0:nyy,2))/ &
	170	((aknew(2)+bknew(2)*psn(0:nxx,0:nyy,1,n,l))- &
	171	(aknew(1)+bknew(1)*psn(0:nxx,0:nyy,1,n,l)))
[db712a8]	172	do kz=2,nz-1
[1a8fbee]	173	pinmconv(0:nxx,0:nyy,kz)= &
	174	(uvzlev(0:nxx,0:nyy,kz+1)-uvzlev(0:nxx,0:nyy,kz-1))/ &
	175	((aknew(kz+1)+bknew(kz+1)*psn(0:nxx,0:nyy,1,n,l))- &
	176	(aknew(kz-1)+bknew(kz-1)*psn(0:nxx,0:nyy,1,n,l)))
[db712a8]	177	end do
[1a8fbee]	178	pinmconv(0:nxx,0:nyy,nz)= &
	179	(uvzlev(0:nxx,0:nyy,nz)-uvzlev(0:nxx,0:nyy,nz-1))/ &
	180	((aknew(nz)+bknew(nz)*psn(0:nxx,0:nyy,1,n,l))- &
	181	(aknew(nz-1)+bknew(nz-1)*psn(0:nxx,0:nyy,1,n,l)))
[db712a8]	182
[1a8fbee]	183	! Levels where u,v,t and q are given
[db712a8]	184	!************************************
	185
[1a8fbee]	186	uun(0:nxx,0:nyy,1,n,l)=uuhn(0:nxx,0:nyy,1,l)
	187	vvn(0:nxx,0:nyy,1,n,l)=vvhn(0:nxx,0:nyy,1,l)
	188	ttn(0:nxx,0:nyy,1,n,l)=tthn(0:nxx,0:nyy,1,n,l)
	189	qvn(0:nxx,0:nyy,1,n,l)=qvhn(0:nxx,0:nyy,1,n,l)
[db712a8]	190	if (readclouds_nest(l)) then
[1a8fbee]	191	clwcn(0:nxx,0:nyy,1,n,l)=clwchn(0:nxx,0:nyy,1,n,l)
	192	if (.not. sumclouds_nest(l)) &
	193	ciwcn(0:nxx,0:nyy,1,n,l)=ciwchn(0:nxx,0:nyy,1,n,l)
[db712a8]	194	end if
[1a8fbee]	195	pvn(0:nxx,0:nyy,1,n,l)=pvhn(0:nxx,0:nyy,1,l)
	196	rhon(0:nxx,0:nyy,1,n,l)=rhohn(0:nxx,0:nyy,1)
[db712a8]	197
[1a8fbee]	198	uun(0:nxx,0:nyy,nz,n,l)=uuhn(0:nxx,0:nyy,nuvz,l)
	199	vvn(0:nxx,0:nyy,nz,n,l)=vvhn(0:nxx,0:nyy,nuvz,l)
	200	ttn(0:nxx,0:nyy,nz,n,l)=tthn(0:nxx,0:nyy,nuvz,n,l)
	201	qvn(0:nxx,0:nyy,nz,n,l)=qvhn(0:nxx,0:nyy,nuvz,n,l)
[db712a8]	202	if (readclouds_nest(l)) then
[1a8fbee]	203	clwcn(0:nxx,0:nyy,nz,n,l)=clwchn(0:nxx,0:nyy,nuvz,n,l)
	204	if (.not. sumclouds_nest(l)) &
	205	ciwcn(0:nxx,0:nyy,nz,n,l)=ciwchn(0:nxx,0:nyy,nuvz,n,l)
[db712a8]	206	end if
[1a8fbee]	207	pvn(0:nxx,0:nyy,nz,n,l)=pvhn(0:nxx,0:nyy,nuvz,l)
	208	rhon(0:nxx,0:nyy,nz,n,l)=rhohn(0:nxx,0:nyy,nuvz)
[db712a8]	209
	210	kmin=2
	211	idx=kmin
[1a8fbee]	212	iz_loop: do iz=2,nz-1
	213
	214	do jy=0,nyy
	215	do ix=0,nxx
	216	if( height(iz).gt.uvzlev(ix,jy,nuvz)) then
	217
[e200b7a]	218	uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
	219	vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
	220	ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
	221	qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
[db712a8]	222	!hg adding the cloud water
	223	if (readclouds_nest(l)) then
	224	clwcn(ix,jy,iz,n,l)=clwcn(ix,jy,nz,n,l)
[1a8fbee]	225	if (.not. sumclouds_nest(l)) &
	226	ciwcn(ix,jy,iz,n,l)=ciwcn(ix,jy,nz,n,l)
[db712a8]	227	end if
	228	!hg
[e200b7a]	229	pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
	230	rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
[1a8fbee]	231
[db712a8]	232	else
[1a8fbee]	233
	234	kz_loop: do kz=idx(ix,jy),nuvz
	235	if (idx(ix,jy) .le. kz .and. height(iz).gt.uvzlev(ix,jy,kz-1) &
	236	.and. height(iz).le.uvzlev(ix,jy,kz)) then
[db712a8]	237	idx(ix,jy)=kz
[1a8fbee]	238	exit kz_loop
[db712a8]	239	endif
[1a8fbee]	240	enddo kz_loop
	241
[e200b7a]	242	endif
[db712a8]	243	enddo
	244	enddo
[1a8fbee]	245
	246	do jy=0,nyy
	247	do ix=0,nxx
[db712a8]	248	if(height(iz).le.uvzlev(ix,jy,nuvz)) then
	249	kz=idx(ix,jy)
	250	dz1=height(iz)-uvzlev(ix,jy,kz-1)
	251	dz2=uvzlev(ix,jy,kz)-height(iz)
[4fbe7a5]	252	dz=dz1+dz2
[db712a8]	253	uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)dz2+uuhn(ix,jy,kz,l)dz1)/dz
	254	vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)dz2+vvhn(ix,jy,kz,l)dz1)/dz
	255	ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2 &
	256	+tthn(ix,jy,kz,n,l)*dz1)/dz
	257	qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2 &
	258	+qvhn(ix,jy,kz,n,l)*dz1)/dz
	259	!hg adding the cloud water
	260	if (readclouds_nest(l)) then
[1a8fbee]	261	clwcn(ix,jy,iz,n,l)=(clwchn(ix,jy,kz-1,n,l)*dz2 + &
	262	clwchn(ix,jy,kz,n,l)*dz1)/dz
	263	if (.not. sumclouds_nest(l)) ciwcn(ix,jy,iz,n,l)= &
	264	(ciwchn(ix,jy,kz-1,n,l)dz2+ciwchn(ix,jy,kz,n,l)dz1)/dz
[db712a8]	265	end if
	266	!hg
	267	pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)dz2+pvhn(ix,jy,kz,l)dz1)/dz
	268	rhon(ix,jy,iz,n,l)=(rhohn(ix,jy,kz-1)dz2+rhohn(ix,jy,kz)dz1)/dz
[1a8fbee]	269
[e200b7a]	270	endif
[db712a8]	271	enddo
	272	enddo
[1a8fbee]	273
	274	enddo iz_loop
[db712a8]	275
	276
	277
	278	! do kz=kmin,nuvz
	279	! if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then
	280	! uun(ix,jy,iz,n,l)=uun(ix,jy,nz,n,l)
	281	! vvn(ix,jy,iz,n,l)=vvn(ix,jy,nz,n,l)
	282	! ttn(ix,jy,iz,n,l)=ttn(ix,jy,nz,n,l)
	283	! qvn(ix,jy,iz,n,l)=qvn(ix,jy,nz,n,l)
	284	! pvn(ix,jy,iz,n,l)=pvn(ix,jy,nz,n,l)
	285	! rhon(ix,jy,iz,n,l)=rhon(ix,jy,nz,n,l)
	286	! goto 30
	287	! endif
	288	! if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
	289	! (height(iz).le.uvwzlev(ix,jy,kz))) then
	290	! dz1=height(iz)-uvwzlev(ix,jy,kz-1)
	291	! dz2=uvwzlev(ix,jy,kz)-height(iz)
	292	! dz=dz1+dz2
	293	! uun(ix,jy,iz,n,l)=(uuhn(ix,jy,kz-1,l)*dz2+ &
	294	! uuhn(ix,jy,kz,l)*dz1)/dz
	295	! vvn(ix,jy,iz,n,l)=(vvhn(ix,jy,kz-1,l)*dz2+ &
	296	! vvhn(ix,jy,kz,l)*dz1)/dz
	297	! ttn(ix,jy,iz,n,l)=(tthn(ix,jy,kz-1,n,l)*dz2+ &
	298	! tthn(ix,jy,kz,n,l)*dz1)/dz
	299	! qvn(ix,jy,iz,n,l)=(qvhn(ix,jy,kz-1,n,l)*dz2+ &
	300	! qvhn(ix,jy,kz,n,l)*dz1)/dz
	301	! pvn(ix,jy,iz,n,l)=(pvhn(ix,jy,kz-1,l)*dz2+ &
	302	! pvhn(ix,jy,kz,l)*dz1)/dz
	303	! rhon(ix,jy,iz,n,l)=(rhohn(kz-1)dz2+rhohn(kz)dz1)/dz
	304	! kmin=kz
	305	! goto 30
	306	! endif
	307	! end do
	308	! 30 continue
	309	! end do
	310
	311
	312	! Levels, where w is given
	313	!*************************
	314
[1a8fbee]	315	wwn(0:nxx,0:nyy,1,n,l)=wwhn(0:nxx,0:nyy,1,l)*pinmconv(0:nxx,0:nyy,1)
	316	wwn(0:nxx,0:nyy,nz,n,l)=wwhn(0:nxx,0:nyy,nwz,l)*pinmconv(0:nxx,0:nyy,nz)
[db712a8]	317	kmin=2
	318	idx=kmin
[1a8fbee]	319	iz_loop2: do iz=2,nz
	320	do jy=0,nyy
	321	do ix=0,nxx
	322	kz_loop2: do kz=idx(ix,jy),nwz
	323	if (idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1) &
	324	.and. height(iz).le.wzlev(ix,jy,kz)) then
[db712a8]	325	idx(ix,jy)=kz
[1a8fbee]	326	exit kz_loop2
[db712a8]	327	endif
[1a8fbee]	328	enddo kz_loop2
[db712a8]	329	enddo
	330	enddo
[1a8fbee]	331	do jy=0,nyy
	332	do ix=0,nxx
[db712a8]	333	kz=idx(ix,jy)
	334	dz1=height(iz)-wzlev(ix,jy,kz-1)
	335	dz2=wzlev(ix,jy,kz)-height(iz)
	336	dz=dz1+dz2
	337	wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)pinmconv(ix,jy,kz-1)dz2 &
	338	+wwhn(ix,jy,kz,l)pinmconv(ix,jy,kz)dz1)/dz
	339	enddo
	340	enddo
[1a8fbee]	341	enddo iz_loop2
[db712a8]	342
	343	! wwn(ix,jy,1,n,l)=wwhn(ix,jy,1,l)*pinmconv(1)
	344	! wwn(ix,jy,nz,n,l)=wwhn(ix,jy,nwz,l)*pinmconv(nz)
	345	! kmin=2
	346	! do iz=2,nz
	347	! do kz=kmin,nwz
	348	! if ((height(iz).gt.wzlev(kz-1)).and. &
	349	! (height(iz).le.wzlev(kz))) then
	350	! dz1=height(iz)-wzlev(kz-1)
	351	! dz2=wzlev(kz)-height(iz)
	352	! dz=dz1+dz2
	353	! wwn(ix,jy,iz,n,l)=(wwhn(ix,jy,kz-1,l)pinmconv(kz-1)dz2 &
	354	! +wwhn(ix,jy,kz,l)pinmconv(kz)dz1)/dz
	355	! kmin=kz
	356	! goto 40
	357	! endif
	358	! end do
	359	! 40 continue
	360	! end do
	361
	362	! Compute density gradients at intermediate levels
	363	!*************************************************
	364
[1a8fbee]	365	drhodzn(0:nxx,0:nyy,1,n,l)= &
	366	( rhon(0:nxx,0:nyy,2,n,l)-rhon(0:nxx,0:nyy,1,n,l) )/(height(2)-height(1))
[db712a8]	367	do kz=2,nz-1
[1a8fbee]	368	drhodzn(0:nxx,0:nyy,kz,n,l)= &
	369	( rhon(0:nxx,0:nyy,kz+1,n,l)-rhon(0:nxx,0:nyy,kz-1,n,l) ) / &
	370	( height(kz+1)-height(kz-1) )
[e200b7a]	371	end do
[1a8fbee]	372	drhodzn(0:nxx,0:nyy,nz,n,l)=drhodzn(0:nxx,0:nyy,nz-1,n,l)
[e200b7a]	373
	374
[db712a8]	375	! drhodzn(ix,jy,1,n,l)=(rhon(ix,jy,2,n,l)-rhon(ix,jy,1,n,l))/ &
	376	! (height(2)-height(1))
	377	! do kz=2,nz-1
	378	! drhodzn(ix,jy,kz,n,l)=(rhon(ix,jy,kz+1,n,l)- &
	379	! rhon(ix,jy,kz-1,n,l))/(height(kz+1)-height(kz-1))
	380	! end do
	381	! drhodzn(ix,jy,nz,n,l)=drhodzn(ix,jy,nz-1,n,l)
[e200b7a]	382
	383
	384
[db712a8]	385	! end do
	386	! end do
[e200b7a]	387
	388
[db712a8]	389	!****************************************************************
	390	! Compute slope of eta levels in windward direction and resulting
	391	! vertical wind correction
	392	!****************************************************************
[e200b7a]	393
[1a8fbee]	394	do jy=1,nyy-1
[db712a8]	395	! cosf=cos((real(jy)dyn(l)+ylat0n(l))pi180)
[1a8fbee]	396	cosf(jy) = 1. / cos( ( real(jy)dyn(l) + ylat0n(l) ) pi180 )
[db712a8]	397	end do
[e200b7a]	398
[db712a8]	399	kmin=2
	400	idx=kmin
[1a8fbee]	401	iz_loop3: do iz=2,nz-1
	402	do jy=1,nyy-1
	403	do ix=1,nxx-1
[db712a8]	404
[1a8fbee]	405	kz_loop3: do kz=idx(ix,jy),nz
	406	if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)) &
	407	.and. (height(iz).le.uvzlev(ix,jy,kz))) then
[db712a8]	408	idx(ix,jy)=kz
[1a8fbee]	409	exit kz_loop3
[db712a8]	410	endif
[1a8fbee]	411	enddo kz_loop3
[db712a8]	412	enddo
	413	enddo
	414
[1a8fbee]	415	do jy=1,nyy-1
	416	do ix=1,nxx-1
[db712a8]	417	kz=idx(ix,jy)
	418	dz1=height(iz)-uvzlev(ix,jy,kz-1)
	419	dz2=uvzlev(ix,jy,kz)-height(iz)
	420	dz=dz1+dz2
	421	ix1=ix-1
	422	jy1=jy-1
	423	ixp=ix+1
	424	jyp=jy+1
	425
	426	dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2.
	427	dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2.
	428	dzdx=(dzdx1dz2+dzdx2dz1)/dz
	429
	430	dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2.
	431	dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2.
	432	dzdy=(dzdy1dz2+dzdy2dz1)/dz
	433
[1a8fbee]	434	wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l) + &
	435	( dzdxuun(ix,jy,iz,n,l)dxconstxresoln(l)cosf(jy) + &
	436	dzdyvvn(ix,jy,iz,n,l)dyconst*yresoln(l) )
[e200b7a]	437
[db712a8]	438	end do
[e200b7a]	439
	440	end do
[1a8fbee]	441	end do iz_loop3
[db712a8]	442
	443
	444	! do jy=1,nyn(l)-2
	445	! do ix=1,nxn(l)-2
	446	! kmin=2
	447	! do iz=2,nz-1
	448
	449	! ui=uun(ix,jy,iz,n,l)dxconstxresoln(l)/cosf(jy)
	450	! vi=vvn(ix,jy,iz,n,l)dyconstyresoln(l)
	451
	452	! do kz=kmin,nz
	453	! if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
	454	! (height(iz).le.uvwzlev(ix,jy,kz))) then
	455	! dz1=height(iz)-uvwzlev(ix,jy,kz-1)
	456	! dz2=uvwzlev(ix,jy,kz)-height(iz)
	457	! dz=dz1+dz2
	458	! kl=kz-1
	459	! klp=kz
	460	! kmin=kz
	461	! goto 47
	462	! endif
	463	! end do
	464
	465	! 47 ix1=ix-1
	466	! jy1=jy-1
	467	! ixp=ix+1
	468	! jyp=jy+1
	469
	470	! dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
	471	! dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
	472	! dzdx=(dzdx1dz2+dzdx2dz1)/dz
	473
	474	! dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
	475	! dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
	476	! dzdy=(dzdy1dz2+dzdy2dz1)/dz
	477
	478	! wwn(ix,jy,iz,n,l)=wwn(ix,jy,iz,n,l)+(dzdxui+dzdyvi)
	479
	480	! end do
	481
	482	! end do
	483	! end do
	484
	485
	486	!write (,) 'initializing nested cloudsn, n:',n
	487	! create a cloud and rainout/washout field, cloudsn occur where rh>80%
	488
	489
[1a8fbee]	490	!******************************************************************************
	491	if (readclouds_nest(l)) then !HG METHOD
	492
	493	! Loops over all grid cells vertically and construct the 3D matrix for clouds
	494	! Cloud top and cloud bottom grid cells are assigned as well as the total column
	495	! cloud water. For precipitating columns, the type and whether it is in or below
[db712a8]	496	! cloud scavenging are assigned with numbers 2-5 (following the old metod).
[1a8fbee]	497	! A distinction is made between lsp and convp though they are treated the equally
	498	! with regard to scavenging. Also, clouds that are not precipitating are defined which
	499	! may be used in the future for cloud processing by non-precipitating-clouds.
	500	!*******************************************************************************
	501
	502	!PS write(,) 'Nested ECMWF fields: using cloud water'
	503	clwn(0:nxx,0:nyy,:,n,l)=0.
	504	! icloud_stats(0:nxx,0:nyy,:,n)=0.
	505	ctwcn(0:nxx,0:nyy,n,l)=0.
	506	cloudsn(0:nxx,0:nyy,:,n,l)=0
[db712a8]	507	! If water/ice are read separately into clwc and ciwc, store sum in clwcn
[1a8fbee]	508	if (.not. sumclouds_nest(l)) then
	509	clwcn(0:nxx,0:nyy,:,n,l) = clwcn(0:nxx,0:nyy,:,n,l) + &
	510	ciwcn(0:nxx,0:nyy,:,n,l)
[db712a8]	511	end if
[1a8fbee]	512	do jy=0,nyy
	513	do ix=0,nxx
[db712a8]	514	lsp=lsprecn(ix,jy,1,n,l)
	515	convp=convprecn(ix,jy,1,n,l)
	516	prec=lsp+convp
	517	tot_cloud_h=0
	518	! Find clouds in the vertical
[1a8fbee]	519	!! Note PS: bad loop order.
[db712a8]	520	do kz=1, nz-1 !go from top to bottom
	521	if (clwcn(ix,jy,kz,n,l).gt.0) then
[1a8fbee]	522	! assuming rho is in kg/m3 and hz in m gives: kg/kg * kg/m3 *m3/kg /m = m2/m3
	523	clwn(ix,jy,kz,n,l)=(clwcn(ix,jy,kz,n,l)rhon(ix,jy,kz,n,l)) &
	524	(height(kz+1)-height(kz))
[db712a8]	525	tot_cloud_h=tot_cloud_h+(height(kz+1)-height(kz))
[341f4b7]	526	ctwcn(ix,jy,n,l) = ctwcn(ix,jy,n,l)+clwn(ix,jy,kz,n,l)
[db712a8]	527	! icloud_stats(ix,jy,4,n)= icloud_stats(ix,jy,4,n)+clw(ix,jy,kz,n) ! Column cloud water [m3/m3]
	528	! icloud_stats(ix,jy,3,n)= min(height(kz+1),height(kz)) ! Cloud BOT height stats [m]
	529	cloudh_min=min(height(kz+1),height(kz))
	530	!ZHG 2015 extra for testing
	531	! clh(ix,jy,kz,n)=height(kz+1)-height(kz)
	532	! icloud_stats(ix,jy,1,n)=icloud_stats(ix,jy,1,n)+(height(kz+1)-height(kz)) ! Cloud total vertical extent [m]
	533	! icloud_stats(ix,jy,2,n)= max(icloud_stats(ix,jy,2,n),height(kz)) ! Cloud TOP height [m]
	534	!ZHG
	535	endif
	536	end do
	537
	538	! If Precipitation. Define removal type in the vertical
[1a8fbee]	539	if (lsp.gt.0.01 .or. convp.gt.0.01) then ! cloud and precipitation
	540	!! Note PS: such hardcoded limits would better be parameters defined in par_mod
[db712a8]	541
	542	do kz=nz,1,-1 !go Bottom up!
[1a8fbee]	543	!! Note PS: bad loop order
	544	if (clwn(ix,jy,kz,n,l) .gt. 0.) then ! is in cloud
	545	cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+height(kz)-height(kz-1)
	546	cloudsn(ix,jy,kz,n,l)=1 ! is a cloud
[db712a8]	547	if (lsp.ge.convp) then
[1a8fbee]	548	cloudsn(ix,jy,kz,n,l)=3 ! lsp in-cloud
[db712a8]	549	else
[1a8fbee]	550	cloudsn(ix,jy,kz,n,l)=2 ! convp in-cloud
	551	endif ! convective or large scale
	552	elseif ( clwn(ix,jy,kz,n,l) .le. 0. .and. &
	553	cloudh_min .ge. height(kz) ) then ! is below cloud
	554	if (lsp .ge. convp) then
	555	cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout
[db712a8]	556	else
[1a8fbee]	557	cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout
	558	endif ! convective or large scale
[db712a8]	559	endif
	560
[1a8fbee]	561	if (height(kz).ge. 19000) then ! set a max height for removal
[db712a8]	562	cloudsn(ix,jy,kz,n,l)=0
[1a8fbee]	563	endif ! clw>0
	564	end do ! kz
[db712a8]	565	endif ! precipitation
	566	end do
	567	end do
[1a8fbee]	568
[db712a8]	569	!********************************************************************
	570	else ! old method:
	571	!********************************************************************
[1a8fbee]	572
	573	!PS write(,) 'Nested fields: using cloud water from Parameterization'
	574	do jy=0,nyy
	575	do ix=0,nxx
[db712a8]	576	rain_cloud_above(ix,jy)=0
	577	lsp=lsprecn(ix,jy,1,n,l)
	578	convp=convprecn(ix,jy,1,n,l)
	579	cloudshn(ix,jy,n,l)=0
	580	do kz_inv=1,nz-1
[1a8fbee]	581	!! Note PS: bad loop order.
[db712a8]	582	kz=nz-kz_inv+1
	583	pressure=rhon(ix,jy,kz,n,l)r_airttn(ix,jy,kz,n,l)
	584	rh=qvn(ix,jy,kz,n,l)/f_qvsat(pressure,ttn(ix,jy,kz,n,l))
	585	cloudsn(ix,jy,kz,n,l)=0
[1a8fbee]	586
	587	if (rh .gt. 0.8) then ! in cloud
	588	!! Note PS: such hardcoded limits would better be parameters defined in par_mod
	589
	590	if (lsp.gt.0.01 .or. convp.gt.0.01) then ! cloud and precipitation
	591	!! Note PS: such hardcoded limits would better be parameters defined in par_mod
[db712a8]	592	rain_cloud_above(ix,jy)=1
	593	cloudshn(ix,jy,n,l)=cloudshn(ix,jy,n,l)+ &
	594	height(kz)-height(kz-1)
	595	if (lsp.ge.convp) then
[e200b7a]	596	cloudsn(ix,jy,kz,n,l)=3 ! lsp dominated rainout
[db712a8]	597	else
[e200b7a]	598	cloudsn(ix,jy,kz,n,l)=2 ! convp dominated rainout
[db712a8]	599	endif
	600	else ! no precipitation
	601	cloudsn(ix,jy,kz,n,l)=1 ! cloud
	602	endif
[1a8fbee]	603
[db712a8]	604	else ! no cloud
[1a8fbee]	605
[db712a8]	606	if (rain_cloud_above(ix,jy).eq.1) then ! scavenging
	607	if (lsp.ge.convp) then
[e200b7a]	608	cloudsn(ix,jy,kz,n,l)=5 ! lsp dominated washout
[db712a8]	609	else
[e200b7a]	610	cloudsn(ix,jy,kz,n,l)=4 ! convp dominated washout
[db712a8]	611	endif
	612	endif
[e200b7a]	613
[1a8fbee]	614	endif
	615	end do ! kz
	616
	617	end do ! ix
	618	end do ! jy
	619
	620	end if! old method:
	621
	622	end do l_loop! end loop over nests
[e200b7a]	623
	624	end subroutine verttransform_nests

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

Download in other formats:

Original Format