Changeset 4d45639 in flexpart.git
- Timestamp:
- Nov 20, 2015, 4:00:21 PM (8 years ago)
- Branches:
- master, 10.4.1_pesei, GFS_025, bugfixes+enhancements, dev, release-10, release-10.4.1, scaling-bug, univie
- Children:
- b255cd0
- Parents:
- adf46ae
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/verttransform.f90
r0f20c31 r4d45639 21 21 22 22 subroutine verttransform(n,uuh,vvh,wwh,pvh) 23 ! i i i i i23 ! i i i i i 24 24 !***************************************************************************** 25 25 ! * … … 67 67 use com_mod 68 68 use cmapf_mod, only: cc2gll 69 ! eso TODO: 70 ! only used for timing of CPU measurement. Remove this (and calls to mpif_mtime below) 71 ! as this routine should not be dependent on MPI 72 ! use mpi_mod 73 ! :TODO 69 74 70 75 implicit none 71 76 72 77 integer :: ix,jy,kz,iz,n,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym 73 integer :: rain_cloud_above ,kz_inv74 ! integer :: icloudtop !PS 75 real :: f_qvsat,pressure76 ! real :: rh,lsp,convp 77 real :: rh,lsp,convp,prec,rhmin,precmin78 real :: rhoh(nuvzmax),pinmconv(nzmax)79 real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi78 integer :: rain_cloud_above(0:nxmax-1,0:nymax-1),kz_inv,idx(0:nxmax-1,0:nymax-1) 79 real :: f_qvsat,pressure,cosf(0:nymax-1) 80 real :: rh,lsp,convp,tim,tim2,rhmin,precmin,prec 81 real :: uvzlev(0:nxmax-1,0:nymax-1,nuvzmax),rhoh(0:nxmax-1,0:nymax-1,nuvzmax) 82 real :: pinmconv(0:nxmax-1,0:nymax-1,nzmax) 83 real :: ew,pint(0:nxmax-1,0:nymax-1),tv(0:nxmax-1,0:nymax-1) 84 real :: tvold(0:nxmax-1,0:nymax-1),pold(0:nxmax-1,0:nymax-1),dz1,dz2,dz,ui,vi 80 85 real :: xlon,ylat,xlonr,dzdx,dzdy 81 real :: dzdx1,dzdx2,dzdy1,dzdy2 ,cosf86 real :: dzdx1,dzdx2,dzdy1,dzdy2 82 87 real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy 83 88 real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax) … … 85 90 real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax) 86 91 real :: wwh(0:nxmax-1,0:nymax-1,nwzmax) 87 real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax) 88 logical lconvectprec 92 real :: wzlev(0:nxmax-1,0:nymax-1,nwzmax) 89 93 real,parameter :: const=r_air/ga 94 ! integer:: ihr, imin, isec, i100th,ihr2, imin2, isec2, i100th2 90 95 parameter (precmin = 0.002) ! minimum prec in mm/h for cloud diagnostics 91 96 … … 103 108 ! If verttransform is called the first time, initialize heights of the * 104 109 ! z levels in meter. The heights are the heights of model levels, where * 105 ! u,v,T and qv are given. * 110 ! u,v,T and qv are given, and of the interfaces, where w is given. So, * 111 ! the vertical resolution in the z system is doubled. As reference point,* 112 ! the lower left corner of the grid is used. * 113 ! Unlike in the eta system, no difference between heights for u,v and * 114 ! heights for w exists. * 106 115 !************************************************************************* 116 117 118 !eso measure CPU time 119 ! call mpif_mtime('verttransform',0) 107 120 108 121 if (init) then … … 110 123 ! Search for a point with high surface pressure (i.e. not above significant topography) 111 124 ! Then, use this point to construct a reference z profile, to be used at all times 112 !***************************************************************************** *********125 !***************************************************************************** 113 126 114 127 do jy=0,nymin1 … … 124 137 125 138 126 tvold =tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &139 tvold(ixm,jym)=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ & 127 140 ps(ixm,jym,1,n)) 128 pold =ps(ixm,jym,1,n)141 pold(ixm,jym)=ps(ixm,jym,1,n) 129 142 height(1)=0. 130 143 … … 133 146 tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n)) 134 147 135 if (abs(tv-tvold).gt.0.2) then 136 height(kz)=height(kz-1)+const*log(pold/pint)* & 137 (tv-tvold)/log(tv/tvold) 148 if (abs(tv(ixm,jym)-tvold(ixm,jym)).gt.0.2) then 149 height(kz)= & 150 height(kz-1)+const*log(pold(ixm,jym)/pint(ixm,jym))* & 151 (tv(ixm,jym)-tvold(ixm,jym))/log(tv(ixm,jym)/tvold(ixm,jym)) 138 152 else 139 height(kz)=height(kz-1)+const*log(pold/pint)*tv 153 height(kz)=height(kz-1)+ & 154 const*log(pold(ixm,jym)/pint(ixm,jym))*tv(ixm,jym) 140 155 endif 141 156 142 tvold =tv143 pold =pint157 tvold(ixm,jym)=tv(ixm,jym) 158 pold(ixm,jym)=pint(ixm,jym) 144 159 end do 145 160 … … 167 182 !************************* 168 183 184 169 185 do jy=0,nymin1 170 186 do ix=0,nxmin1 171 tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ps(ix,jy,1,n)) 172 pold=ps(ix,jy,1,n) 173 uvwzlev(ix,jy,1)=0. 174 wzlev(1)=0. 175 rhoh(1)=pold/(r_air*tvold) 187 tvold(ix,jy)=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ & 188 ps(ix,jy,1,n)) 189 enddo 190 enddo 191 pold=ps(:,:,1,n) 192 uvzlev(:,:,1)=0. 193 wzlev(:,:,1)=0. 194 rhoh(:,:,1)=pold/(r_air*tvold) 176 195 177 196 … … 179 198 !****************************** 180 199 181 do kz=2,nuvz 182 pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n) 183 tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n)) 184 rhoh(kz)=pint/(r_air*tv) 185 186 if (abs(tv-tvold).gt.0.2) then 187 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)* & 188 (tv-tvold)/log(tv/tvold) 189 else 190 uvwzlev(ix,jy,kz)=uvwzlev(ix,jy,kz-1)+const*log(pold/pint)*tv 191 endif 192 193 tvold=tv 194 pold=pint 195 end do 196 197 198 do kz=2,nwz-1 199 wzlev(kz)=(uvwzlev(ix,jy,kz+1)+uvwzlev(ix,jy,kz))/2. 200 end do 201 wzlev(nwz)=wzlev(nwz-1)+uvwzlev(ix,jy,nuvz)-uvwzlev(ix,jy,nuvz-1) 200 do kz=2,nuvz 201 pint=akz(kz)+bkz(kz)*ps(:,:,1,n) 202 tv=tth(:,:,kz,n)*(1.+0.608*qvh(:,:,kz,n)) 203 rhoh(:,:,kz)=pint(:,:)/(r_air*tv) 204 205 where (abs(tv-tvold).gt.0.2) 206 uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)* & 207 (tv-tvold)/log(tv/tvold) 208 elsewhere 209 uvzlev(:,:,kz)=uvzlev(:,:,kz-1)+const*log(pold/pint)*tv 210 endwhere 211 212 tvold=tv 213 pold=pint 214 end do 215 216 217 do kz=2,nwz-1 218 wzlev(:,:,kz)=(uvzlev(:,:,kz+1)+uvzlev(:,:,kz))/2. 219 end do 220 wzlev(:,:,nwz)=wzlev(:,:,nwz-1)+ & 221 uvzlev(:,:,nuvz)-uvzlev(:,:,nuvz-1) 202 222 203 223 ! pinmconv=(h2-h1)/(p2-p1) 204 224 205 pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &206 ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &207 (aknew(1)+bknew(1)*ps(ix,jy,1,n)))208 209 pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &210 ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &211 (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))212 213 pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &214 ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &215 (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))225 pinmconv(:,:,1)=(uvzlev(:,:,2))/ & 226 ((aknew(2)+bknew(2)*ps(:,:,1,n))- & 227 (aknew(1)+bknew(1)*ps(:,:,1,n))) 228 do kz=2,nz-1 229 pinmconv(:,:,kz)=(uvzlev(:,:,kz+1)-uvzlev(:,:,kz-1))/ & 230 ((aknew(kz+1)+bknew(kz+1)*ps(:,:,1,n))- & 231 (aknew(kz-1)+bknew(kz-1)*ps(:,:,1,n))) 232 end do 233 pinmconv(:,:,nz)=(uvzlev(:,:,nz)-uvzlev(:,:,nz-1))/ & 234 ((aknew(nz)+bknew(nz)*ps(:,:,1,n))- & 235 (aknew(nz-1)+bknew(nz-1)*ps(:,:,1,n))) 216 236 217 237 ! Levels, where u,v,t and q are given 218 238 !************************************ 219 239 220 uu(ix,jy,1,n)=uuh(ix,jy,1) 221 vv(ix,jy,1,n)=vvh(ix,jy,1) 222 tt(ix,jy,1,n)=tth(ix,jy,1,n) 223 qv(ix,jy,1,n)=qvh(ix,jy,1,n) 240 241 uu(:,:,1,n)=uuh(:,:,1) 242 vv(:,:,1,n)=vvh(:,:,1) 243 tt(:,:,1,n)=tth(:,:,1,n) 244 qv(:,:,1,n)=qvh(:,:,1,n) 224 245 !hg adding the cloud water 225 clwc(ix,jy,1,n)=clwch(ix,jy,1,n)226 ciwc(ix,jy,1,n)=ciwch(ix,jy,1,n)246 clwc(:,:,1,n)=clwch(:,:,1,n) 247 ciwc(:,:,1,n)=ciwch(:,:,1,n) 227 248 !hg 228 pv(ix,jy,1,n)=pvh(ix,jy,1)229 rho(ix,jy,1,n)=rhoh(1)230 uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)231 vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)232 tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)233 qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)249 pv(:,:,1,n)=pvh(:,:,1) 250 rho(:,:,1,n)=rhoh(:,:,1) 251 uu(:,:,nz,n)=uuh(:,:,nuvz) 252 vv(:,:,nz,n)=vvh(:,:,nuvz) 253 tt(:,:,nz,n)=tth(:,:,nuvz,n) 254 qv(:,:,nz,n)=qvh(:,:,nuvz,n) 234 255 235 256 !hg adding the cloud water 236 clwc(ix,jy,nz,n)=clwch(ix,jy,nuvz,n)237 ciwc(ix,jy,nz,n)=ciwch(ix,jy,nuvz,n)257 clwc(:,:,nz,n)=clwch(:,:,nuvz,n) 258 ciwc(:,:,nz,n)=ciwch(:,:,nuvz,n) 238 259 !hg 239 pv(ix,jy,nz,n)=pvh(ix,jy,nuvz) 240 rho(ix,jy,nz,n)=rhoh(nuvz) 241 kmin=2 242 do iz=2,nz-1 243 do kz=kmin,nuvz 244 if(height(iz).gt.uvwzlev(ix,jy,nuvz)) then 245 uu(ix,jy,iz,n)=uu(ix,jy,nz,n) 246 vv(ix,jy,iz,n)=vv(ix,jy,nz,n) 247 tt(ix,jy,iz,n)=tt(ix,jy,nz,n) 248 qv(ix,jy,iz,n)=qv(ix,jy,nz,n) 260 pv(:,:,nz,n)=pvh(:,:,nuvz) 261 rho(:,:,nz,n)=rhoh(:,:,nuvz) 262 263 264 kmin=2 265 idx=kmin 266 do iz=2,nz-1 267 do jy=0,nymin1 268 do ix=0,nxmin1 269 if(height(iz).gt.uvzlev(ix,jy,nuvz)) then 270 uu(ix,jy,iz,n)=uu(ix,jy,nz,n) 271 vv(ix,jy,iz,n)=vv(ix,jy,nz,n) 272 tt(ix,jy,iz,n)=tt(ix,jy,nz,n) 273 qv(ix,jy,iz,n)=qv(ix,jy,nz,n) 249 274 !hg adding the cloud water 250 251 275 clwc(ix,jy,iz,n)=clwc(ix,jy,nz,n) 276 ciwc(ix,jy,iz,n)=ciwc(ix,jy,nz,n) 252 277 !hg 253 pv(ix,jy,iz,n)=pv(ix,jy,nz,n) 254 rho(ix,jy,iz,n)=rho(ix,jy,nz,n) 255 goto 30 256 endif 257 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 258 (height(iz).le.uvwzlev(ix,jy,kz))) then 259 dz1=height(iz)-uvwzlev(ix,jy,kz-1) 260 dz2=uvwzlev(ix,jy,kz)-height(iz) 261 dz=dz1+dz2 262 uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz 263 vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz 264 tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2+tth(ix,jy,kz,n)*dz1)/dz 265 qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2+qvh(ix,jy,kz,n)*dz1)/dz 278 pv(ix,jy,iz,n)=pv(ix,jy,nz,n) 279 rho(ix,jy,iz,n)=rho(ix,jy,nz,n) 280 else 281 innuvz: do kz=idx(ix,jy),nuvz 282 if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. & 283 (height(iz).le.uvzlev(ix,jy,kz))) then 284 idx(ix,jy)=kz 285 exit innuvz 286 endif 287 enddo innuvz 288 endif 289 enddo 290 enddo 291 do jy=0,nymin1 292 do ix=0,nxmin1 293 if(height(iz).le.uvzlev(ix,jy,nuvz)) then 294 kz=idx(ix,jy) 295 dz1=height(iz)-uvzlev(ix,jy,kz-1) 296 dz2=uvzlev(ix,jy,kz)-height(iz) 297 dz=dz1+dz2 298 uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz 299 vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz 300 tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 & 301 +tth(ix,jy,kz,n)*dz1)/dz 302 qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 & 303 +qvh(ix,jy,kz,n)*dz1)/dz 266 304 !hg adding the cloud water 267 268 305 clwc(ix,jy,iz,n)=(clwch(ix,jy,kz-1,n)*dz2+clwch(ix,jy,kz,n)*dz1)/dz 306 ciwc(ix,jy,iz,n)=(ciwch(ix,jy,kz-1,n)*dz2+ciwch(ix,jy,kz,n)*dz1)/dz 269 307 !hg 270 271 pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz 272 rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz 273 kmin=kz 274 goto 30 275 endif 276 end do 277 30 continue 278 end do 308 pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz 309 rho(ix,jy,iz,n)=(rhoh(ix,jy,kz-1)*dz2+rhoh(ix,jy,kz)*dz1)/dz 310 endif 311 enddo 312 enddo 313 enddo 279 314 280 315 … … 282 317 !************************* 283 318 284 ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1) 285 ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz) 286 kmin=2 287 do iz=2,nz 288 do kz=kmin,nwz 289 if ((height(iz).gt.wzlev(kz-1)).and. & 290 (height(iz).le.wzlev(kz))) then 291 dz1=height(iz)-wzlev(kz-1) 292 dz2=wzlev(kz)-height(iz) 293 dz=dz1+dz2 294 ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 & 295 +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz 296 kmin=kz 297 goto 40 319 ww(:,:,1,n)=wwh(:,:,1)*pinmconv(:,:,1) 320 ww(:,:,nz,n)=wwh(:,:,nwz)*pinmconv(:,:,nz) 321 kmin=2 322 idx=kmin 323 do iz=2,nz 324 do jy=0,nymin1 325 do ix=0,nxmin1 326 inn: do kz=idx(ix,jy),nwz 327 if(idx(ix,jy) .le. kz .and. height(iz).gt.wzlev(ix,jy,kz-1).and. & 328 height(iz).le.wzlev(ix,jy,kz)) then 329 idx(ix,jy)=kz 330 exit inn 298 331 endif 299 end do 300 40 continue 301 end do 332 enddo inn 333 enddo 334 enddo 335 do jy=0,nymin1 336 do ix=0,nxmin1 337 kz=idx(ix,jy) 338 dz1=height(iz)-wzlev(ix,jy,kz-1) 339 dz2=wzlev(ix,jy,kz)-height(iz) 340 dz=dz1+dz2 341 ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(ix,jy,kz-1)*dz2 & 342 +wwh(ix,jy,kz)*pinmconv(ix,jy,kz)*dz1)/dz 343 enddo 344 enddo 345 enddo 302 346 303 347 ! Compute density gradients at intermediate levels 304 348 !************************************************* 305 349 306 drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ & 307 (height(2)-height(1)) 308 do kz=2,nz-1 309 drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ & 310 (height(kz+1)-height(kz-1)) 311 end do 312 drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n) 313 314 end do 350 drhodz(:,:,1,n)=(rho(:,:,2,n)-rho(:,:,1,n))/ & 351 (height(2)-height(1)) 352 do kz=2,nz-1 353 drhodz(:,:,kz,n)=(rho(:,:,kz+1,n)-rho(:,:,kz-1,n))/ & 354 (height(kz+1)-height(kz-1)) 315 355 end do 356 drhodz(:,:,nz,n)=drhodz(:,:,nz-1,n) 357 358 ! end do 359 ! end do 316 360 317 361 … … 322 366 323 367 do jy=1,ny-2 324 cosf=cos((real(jy)*dy+ylat0)*pi180) 325 do ix=1,nx-2 326 327 kmin=2 328 do iz=2,nz-1 329 330 ui=uu(ix,jy,iz,n)*dxconst/cosf 331 vi=vv(ix,jy,iz,n)*dyconst 332 333 do kz=kmin,nz 334 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. & 335 (height(iz).le.uvwzlev(ix,jy,kz))) then 336 dz1=height(iz)-uvwzlev(ix,jy,kz-1) 337 dz2=uvwzlev(ix,jy,kz)-height(iz) 338 dz=dz1+dz2 339 kl=kz-1 340 klp=kz 341 kmin=kz 342 goto 47 368 cosf(jy)=1./cos((real(jy)*dy+ylat0)*pi180) 369 enddo 370 371 kmin=2 372 idx=kmin 373 do iz=2,nz-1 374 do jy=1,ny-2 375 do ix=1,nx-2 376 377 inneta: do kz=idx(ix,jy),nz 378 if (idx(ix,jy) .le. kz .and. (height(iz).gt.uvzlev(ix,jy,kz-1)).and. & 379 (height(iz).le.uvzlev(ix,jy,kz))) then 380 idx(ix,jy)=kz 381 exit inneta 343 382 endif 344 end do 345 346 47 ix1=ix-1 383 enddo inneta 384 enddo 385 enddo 386 387 do jy=1,ny-2 388 do ix=1,nx-2 389 kz=idx(ix,jy) 390 dz1=height(iz)-uvzlev(ix,jy,kz-1) 391 dz2=uvzlev(ix,jy,kz)-height(iz) 392 dz=dz1+dz2 393 ix1=ix-1 347 394 jy1=jy-1 348 395 ixp=ix+1 349 396 jyp=jy+1 350 397 351 dzdx1=(uv wzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.352 dzdx2=(uv wzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.398 dzdx1=(uvzlev(ixp,jy,kz-1)-uvzlev(ix1,jy,kz-1))/2. 399 dzdx2=(uvzlev(ixp,jy,kz)-uvzlev(ix1,jy,kz))/2. 353 400 dzdx=(dzdx1*dz2+dzdx2*dz1)/dz 354 401 355 dzdy1=(uv wzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.356 dzdy2=(uv wzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.402 dzdy1=(uvzlev(ix,jyp,kz-1)-uvzlev(ix,jy1,kz-1))/2. 403 dzdy2=(uvzlev(ix,jyp,kz)-uvzlev(ix,jy1,kz))/2. 357 404 dzdy=(dzdy1*dz2+dzdy2*dz1)/dz 358 405 359 ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*u i+dzdy*vi)406 ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*uu(ix,jy,iz,n)*dxconst*cosf(jy)+dzdy*vv(ix,jy,iz,n)*dyconst) 360 407 361 408 end do … … 363 410 end do 364 411 end do 365 366 412 367 413 ! If north pole is in the domain, calculate wind velocities in polar … … 370 416 371 417 if (nglobal) then 372 do jy=int(switchnorthg)-2,nymin1373 ylat=ylat0+real(jy)*dy374 do ix=0,nxmin1375 xlon=xlon0+real(ix)*dx376 do iz=1,nz418 do iz=1,nz 419 do jy=int(switchnorthg)-2,nymin1 420 ylat=ylat0+real(jy)*dy 421 do ix=0,nxmin1 422 xlon=xlon0+real(ix)*dx 377 423 call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), & 378 vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n)) 424 vv(ix,jy,iz,n),uupol(ix,jy,iz,n), & 425 vvpol(ix,jy,iz,n)) 379 426 end do 380 427 end do … … 390 437 xlon=xlon0+real(nx/2-1)*dx 391 438 xlonr=xlon*pi/180. 392 ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+vv(nx/2-1,nymin1,iz,n)**2) 439 ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ & 440 vv(nx/2-1,nymin1,iz,n)**2) 393 441 if (vv(nx/2-1,nymin1,iz,n).lt.0.) then 394 ddpol=atan(uu(nx/2-1,nymin1,iz,n)/vv(nx/2-1,nymin1,iz,n))-xlonr 442 ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ & 443 vv(nx/2-1,nymin1,iz,n))-xlonr 395 444 else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then 396 445 ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ & … … 408 457 uuaux=-ffpol*sin(xlonr+ddpol) 409 458 vvaux=-ffpol*cos(xlonr+ddpol) 410 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux) 459 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, & 460 vvpolaux) 461 411 462 jy=nymin1 412 463 do ix=0,nxmin1 … … 441 492 442 493 if (sglobal) then 443 do jy=0,int(switchsouthg)+3444 ylat=ylat0+real(jy)*dy445 do ix=0,nxmin1446 xlon=xlon0+real(ix)*dx447 do iz=1,nz494 do iz=1,nz 495 do jy=0,int(switchsouthg)+3 496 ylat=ylat0+real(jy)*dy 497 do ix=0,nxmin1 498 xlon=xlon0+real(ix)*dx 448 499 call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), & 449 vv(ix,jy,iz,n),uupol(ix,jy,iz,n),vvpol(ix,jy,iz,n)) 500 vv(ix,jy,iz,n),uupol(ix,jy,iz,n), & 501 vvpol(ix,jy,iz,n)) 450 502 end do 451 503 end do … … 460 512 xlon=xlon0+real(nx/2-1)*dx 461 513 xlonr=xlon*pi/180. 462 ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+vv(nx/2-1,0,iz,n)**2) 514 ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ & 515 vv(nx/2-1,0,iz,n)**2) 463 516 if (vv(nx/2-1,0,iz,n).lt.0.) then 464 ddpol=atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr 517 ddpol=atan(uu(nx/2-1,0,iz,n)/ & 518 vv(nx/2-1,0,iz,n))+xlonr 465 519 else if (vv(nx/2-1,0,iz,n).gt.0.) then 466 ddpol=pi+atan(uu(nx/2-1,0,iz,n)/vv(nx/2-1,0,iz,n))+xlonr 520 ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ & 521 vv(nx/2-1,0,iz,n))+xlonr 467 522 else 468 523 ddpol=pi/2-xlonr … … 477 532 uuaux=+ffpol*sin(xlonr-ddpol) 478 533 vvaux=-ffpol*cos(xlonr-ddpol) 479 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux,vvpolaux) 534 call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, & 535 vvpolaux) 480 536 481 537 jy=0 … … 505 561 506 562 507 !write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz508 ! create a cloud and rainout/washout field, clouds occur where rh>80%509 ! total cloudheight is stored at level 0510 511 563 if (readclouds) write(*,*) 'using cloud water from ECMWF' 512 564 ! if (.not.readclouds) write(*,*) 'using cloud water from parameterization' … … 518 570 do ix=0,nxmin1 519 571 ! OLD METHOD 520 rain_cloud_above =0572 rain_cloud_above(ix,jy)=0 521 573 lsp=lsprec(ix,jy,1,n) 522 574 convp=convprec(ix,jy,1,n) … … 529 581 if (rh.gt.0.8) then ! in cloud 530 582 if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation 531 rain_cloud_above =1583 rain_cloud_above(ix,jy)=1 532 584 cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ & 533 585 height(kz)-height(kz-1) … … 541 593 endif 542 594 else ! no cloud 543 if (rain_cloud_above .eq.1) then ! scavenging595 if (rain_cloud_above(ix,jy).eq.1) then ! scavenging 544 596 if (lsp.ge.convp) then 545 597 clouds(ix,jy,kz,n)=5 ! lsp dominated washout … … 553 605 554 606 ! PS 2012 555 ! 556 ! 607 ! lsp=lsprec(ix,jy,1,n) 608 ! convp=convprec(ix,jy,1,n) 557 609 ! prec=lsp+convp 558 610 ! if (lsp.gt.convp) then ! prectype='lsp' … … 560 612 ! else ! prectype='cp ' 561 613 ! lconvectprec = .true. 562 ! endif614 ! endif 563 615 !HG METHOD 564 616 !readclouds =.true. … … 572 624 ! do kz=1, nz 573 625 ! !clw & ciw are given in kg/kg 574 ! cloud_water=clwc(ix,jy,kz,n)+ciwc(ix,jy,kz,n) 575 ! if (cloud_water .gt. 0) then 576 ! cloud_min=min(nint(height(kz)),cloud_min) !hg needs reset each grid 577 ! cloud_max=max(nint(height(kz)),cloud_max) !hg needs reset each grid 578 ! cloud_col_wat=cloud_col_wat+cloud_water !hg needs reset each grid 579 ! endif 580 ! cloud_ver=max(0,cloud_max-cloud_min) 581 582 ! if (clwc(ix,jy,kz,n).gt.0 .or. ciwc(ix,jy,kz,n).gt.0) & 583 ! !write(*,*) 'WATER',clwc(ix,jy,kz,n), 'ICE',ciwc(ix,jy,kz,n),'RH',rh,'KZ',kz & 584 ! write(*,*) 'WATER',cloud_water,'RH',rh,'PREC',prec,'HEIGHT',height(kz) & 585 ! enddo 586 ! if ( cloud_min .ne. 99999 .and. cloud_max .ne. -1) write(*,*) 'CB', cloud_min, ' CT',cloud_max 587 ! if (prec .gt. 0) write(*,*) 'PREC',prec,'Cloud Bot',cloud_min,'Cloud Top',cloud_max, 'Cloud Vert. ext',cloud_ver & 588 ! ,'COLUMN cloud water',cloud_col_wat,'Conevctive' ,lconvectprec 589 ! icloudbot(ix,jy,n)=cloud_min 590 ! icloudthck(ix,jy,n)=cloud_ver 591 ! rcw(ix,jy)=cloud_col_wat 592 ! rpc(ix,jy)=prec 626 ! cloud_water=clwc(ix,jy,kz,n)+ciwc(ix,jy,kz,n) 627 ! if (cloud_water .gt. 0) then 628 ! cloud_min=min(nint(height(kz)),cloud_min) !hg needs reset each grid 629 ! cloud_max=max(nint(height(kz)),cloud_max) !hg needs reset each grid 630 ! cloud_col_wat=cloud_col_wat+cloud_water !hg needs reset each grid 631 ! endif 632 ! cloud_ver=max(0,cloud_max-cloud_min) 633 634 ! icloudbot(ix,jy,n)=cloud_min 635 ! icloudthck(ix,jy,n)=cloud_ver 636 ! rcw(ix,jy)=cloud_col_wat 637 ! rpc(ix,jy)=prec 593 638 !write(*,*) 'Using clouds from ECMWF' !hg END Henrik Code 594 639 !END HG METHOD … … 613 658 ! icloudtop=nint(height(kz)) ! use int to save memory 614 659 ! endif 615 ! enddo660 ! end do 616 661 !PS try to get a cloud thicker than 50 m 617 662 !PS if there is at least .01 mm/h - changed to 0.002 and put into … … 622 667 ! rhmin = rhmin - 0.05 623 668 ! if (rhmin .ge. 0.30) goto 98 ! give up for <= 25% rel.hum. 624 ! endif669 ! end if 625 670 626 671 !PS is based on looking at a limited set of comparison data … … 656 701 end do 657 702 658 end subroutine verttransform 703 !do 102 kz=1,nuvz 704 !write(an,'(i02)') kz+10 705 !write(*,*) nuvz,nymin1,nxmin1,'--',an,'--' 706 !open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted') 707 !do 101 jy=0,nymin1 708 ! write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1) 709 !101 continue 710 ! close(4) 711 !102 continue 712 713 ! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted') 714 ! do 103 jy=0,nymin1 715 ! write (4,*) 716 !+ (height(kz),kz=1,nuvz) 717 !103 continue 718 ! close(4) 719 720 !open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted') 721 ! do 104 jy=0,nymin1 722 ! write (4,*) 723 !+ (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1) 724 !104 continue 725 ! close(4) 726 727 !eso measure CPU time 728 ! call mpif_mtime('verttransform',1) 729 730 !eso print out the same measure as in Leo's routine 731 ! write(*,*) 'verttransform: ', & 732 ! sum(tt(:,:,:,n)*tt(:,:,:,n)), & 733 ! sum(uu(:,:,:,n)*uu(:,:,:,n)),sum(vv(:,:,:,n)*vv(:,:,:,n)),& 734 ! sum(qv(:,:,:,n)*qv(:,:,:,n)),sum(pv(:,:,:,n)*pv(:,:,:,n)),& 735 ! sum(rho(:,:,:,n)*rho(:,:,:,n)),sum(drhodz(:,:,:,n)*drhodz(:,:,:,n)),& 736 ! sum(ww(:,:,:,n)*ww(:,:,:,n)), & 737 ! sum(clouds(:,:,:,n)), sum(cloudsh(:,:,n)),sum(idx),sum(pinmconv) 738 end subroutine verttransform 739
Note: See TracChangeset
for help on using the changeset viewer.