69 real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
70 real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
71 real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
72 real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
106 integer :: rain_cloud_above,kz_inv
109 real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
110 real ::
ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
111 real :: xlon,ylat,xlonr,dzdx,dzdy
112 real :: dzdx1,dzdx2,dzdy1,dzdy2
113 real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
114 real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
120 integer :: ix,jy,kz,iz,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
125 real,
parameter :: const=r_air/ga
161 logical :: init = .true.
187 if (ps(ix,jy,1,n).gt.100000.)
then
197 tvold=tt2(ixm,jym,1,n)*(1.+0.378*
ew(td2(ixm,jym,1,n))/ &
203 pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
204 tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
216 if (abs(tv-tvold).gt.0.2)
then
218 height(kz-1)+const*log(pold/pint)* &
219 (tv-tvold)/log(tv/tvold)
221 height(kz)=height(kz-1)+ &
222 const*log(pold/pint)*tv
254 if (height(kz).gt.hmixmax)
then
274 tvold=tt2(ix,jy,1,n)*(1.+0.378*
ew(td2(ix,jy,1,n))/ &
279 rhoh(1)=pold/(r_air*tvold)
286 pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
287 tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
288 rhoh(kz)=pint/(r_air*tv)
290 if (abs(tv-tvold).gt.0.2)
then
291 uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
292 (tv-tvold)/log(tv/tvold)
294 uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
303 wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
305 wzlev(nwz)=wzlev(nwz-1)+ &
306 uvzlev(nuvz)-uvzlev(nuvz-1)
310 uvwzlev(ix,jy,kz)=uvzlev(kz)
323 pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
324 ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
325 (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
327 pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
328 ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
329 (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
331 pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
332 ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
333 (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
338 uu(ix,jy,1,n)=uuh(ix,jy,1)
339 vv(ix,jy,1,n)=vvh(ix,jy,1)
340 tt(ix,jy,1,n)=tth(ix,jy,1,n)
341 qv(ix,jy,1,n)=qvh(ix,jy,1,n)
342 pv(ix,jy,1,n)=pvh(ix,jy,1)
343 rho(ix,jy,1,n)=rhoh(1)
344 uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
345 vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
346 tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
347 qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
348 pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
349 rho(ix,jy,nz,n)=rhoh(nuvz)
353 if(height(iz).gt.uvzlev(nuvz))
then
354 uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
355 vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
356 tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
357 qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
358 pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
359 rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
362 if ((height(iz).gt.uvzlev(kz-1)).and. &
363 (height(iz).le.uvzlev(kz)))
then
364 dz1=height(iz)-uvzlev(kz-1)
365 dz2=uvzlev(kz)-height(iz)
367 uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
368 vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
369 tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
370 +tth(ix,jy,kz,n)*dz1)/dz
371 qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
372 +qvh(ix,jy,kz,n)*dz1)/dz
373 pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
374 rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
386 ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
387 ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
391 if ((height(iz).gt.wzlev(kz-1)).and. &
392 (height(iz).le.wzlev(kz)))
then
393 dz1=height(iz)-wzlev(kz-1)
394 dz2=wzlev(kz)-height(iz)
396 ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
397 +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
408 drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
409 (height(2)-height(1))
411 drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
412 (height(kz+1)-height(kz-1))
414 drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
431 ui=uu(ix,jy,iz,n)*dxconst/cos((
real(jy)*dy+ylat0)*pi180)
432 vi=vv(ix,jy,iz,n)*dyconst
435 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
436 (height(iz).le.uvwzlev(ix,jy,kz)))
then
437 dz1=height(iz)-uvwzlev(ix,jy,kz-1)
438 dz2=uvwzlev(ix,jy,kz)-height(iz)
452 dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
453 dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
454 dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
456 dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
457 dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
458 dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
460 ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)
473 do jy=int(switchnorthg)-2,nymin1
474 ylat=ylat0+
real(jy)*dy
476 xlon=xlon0+
real(ix)*dx
478 call
cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
479 vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
492 xlon=xlon0+
real(nx/2-1)*dx
494 ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
495 vv(nx/2-1,nymin1,iz,n)**2)
496 if (vv(nx/2-1,nymin1,iz,n).lt.0.)
then
497 ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
498 vv(nx/2-1,nymin1,iz,n))-xlonr
499 else if (vv(nx/2-1,nymin1,iz,n).gt.0.)
then
500 ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
501 vv(nx/2-1,nymin1,iz,n))-xlonr
505 if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
506 if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
512 uuaux=-ffpol*sin(xlonr+ddpol)
513 vvaux=-ffpol*cos(xlonr+ddpol)
514 call
cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
519 uupol(ix,jy,iz,n)=uupolaux
520 vvpol(ix,jy,iz,n)=vvpolaux
532 wdummy=wdummy+ww(ix,jy,iz,n)
534 wdummy=wdummy/
real(nx)
537 ww(ix,jy,iz,n)=wdummy
549 do jy=0,int(switchsouthg)+3
550 ylat=ylat0+
real(jy)*dy
552 xlon=xlon0+
real(ix)*dx
554 call
cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
555 vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
567 xlon=xlon0+
real(nx/2-1)*dx
569 ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
570 vv(nx/2-1,0,iz,n)**2)
571 if (vv(nx/2-1,0,iz,n).lt.0.)
then
572 ddpol=atan(uu(nx/2-1,0,iz,n)/ &
573 vv(nx/2-1,0,iz,n))+xlonr
574 else if (vv(nx/2-1,0,iz,n).gt.0.)
then
575 ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
576 vv(nx/2-1,0,iz,n))+xlonr
580 if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
581 if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
587 uuaux=+ffpol*sin(xlonr-ddpol)
588 vvaux=-ffpol*cos(xlonr-ddpol)
589 call
cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
594 uupol(ix,jy,iz,n)=uupolaux
595 vvpol(ix,jy,iz,n)=vvpolaux
607 wdummy=wdummy+ww(ix,jy,iz,n)
609 wdummy=wdummy/
real(nx)
612 ww(ix,jy,iz,n)=wdummy
624 lsp=lsprec(ix,jy,1,n)
625 convp=convprec(ix,jy,1,n)
629 pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
630 rh=qv(ix,jy,kz,n)/
f_qvsat(pressure,tt(ix,jy,kz,n))
633 if ((lsp.gt.0.01).or.(convp.gt.0.01))
then
635 cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
636 height(kz)-height(kz-1)
637 if (lsp.ge.convp)
then
646 if (rain_cloud_above.eq.1)
then
647 if (lsp.ge.convp)
then
subroutine, public cc2gll(strcmp, xlat, xlong, ue, vn, ug, vg)
real function f_qvsat(p, t)