66 real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
67 real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
68 real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
69 real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
103 integer :: rain_cloud_above,kz_inv
106 real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
107 real ::
ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
108 real :: xlon,ylat,xlonr,dzdx,dzdy
109 real :: dzdx1,dzdx2,dzdy1,dzdy2
110 real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
111 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.
182 if (ps(ix,jy,1,n).gt.100000.)
then
192 tvold=tt2(ixm,jym,1,n)*(1.+0.378*
ew(td2(ixm,jym,1,n))/ &
198 pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
199 tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
211 if (abs(tv-tvold).gt.0.2)
then
213 height(kz-1)+const*log(pold/pint)* &
214 (tv-tvold)/log(tv/tvold)
216 height(kz)=height(kz-1)+ &
217 const*log(pold/pint)*tv
249 if (height(kz).gt.hmixmax)
then
273 if (ps(ix,jy,1,n).lt.akz(i)) llev=i
276 if (llev.gt.nuvz-2) llev = nuvz-2
285 tvold=tth(ix,jy,llev,n)*(1.+0.608*qvh(ix,jy,llev,n))
289 uvwzlev(ix,jy,llev)=0.
290 rhoh(llev)=pold/(r_air*tvold)
293 pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
294 tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
295 rhoh(kz)=pint/(r_air*tv)
297 if (abs(tv-tvold).gt.0.2)
then
298 uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
299 (tv-tvold)/log(tv/tvold)
301 uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
304 uvwzlev(ix,jy,kz)=uvzlev(kz)
321 pinmconv(llev)=(uvwzlev(ix,jy,llev+1)-uvwzlev(ix,jy,llev))/ &
322 ((aknew(llev+1)+bknew(llev+1)*ps(ix,jy,1,n))- &
323 (aknew(llev)+bknew(llev)*ps(ix,jy,1,n)))
325 pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
326 ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
327 (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
329 pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
330 ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
331 (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
337 uu(ix,jy,1,n)=uuh(ix,jy,llev)
338 vv(ix,jy,1,n)=vvh(ix,jy,llev)
339 tt(ix,jy,1,n)=tth(ix,jy,llev,n)
340 qv(ix,jy,1,n)=qvh(ix,jy,llev,n)
341 pv(ix,jy,1,n)=pvh(ix,jy,llev)
342 rho(ix,jy,1,n)=rhoh(llev)
343 pplev(ix,jy,1,n)=akz(llev)
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)
350 pplev(ix,jy,nz,n)=akz(nuvz)
354 if(height(iz).gt.uvzlev(nuvz))
then
355 uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
356 vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
357 tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
358 qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
359 pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
360 rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
361 pplev(ix,jy,iz,n)=pplev(ix,jy,nz,n)
364 if ((height(iz).gt.uvzlev(kz-1)).and. &
365 (height(iz).le.uvzlev(kz)))
then
366 dz1=height(iz)-uvzlev(kz-1)
367 dz2=uvzlev(kz)-height(iz)
369 uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
370 vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
371 tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
372 +tth(ix,jy,kz,n)*dz1)/dz
373 qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
374 +qvh(ix,jy,kz,n)*dz1)/dz
375 pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
376 rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
377 pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz
387 ww(ix,jy,1,n)=wwh(ix,jy,llev)*pinmconv(llev)
388 ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
392 if ((height(iz).gt.wzlev(kz-1)).and. &
393 (height(iz).le.wzlev(kz)))
then
394 dz1=height(iz)-wzlev(kz-1)
395 dz2=wzlev(kz)-height(iz)
397 ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
398 +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 if (ps(ix,jy,1,n).lt.akz(i)) llev=i
434 if (llev.gt.nuvz-2) llev = nuvz-2
442 ui=uu(ix,jy,iz,n)*dxconst/cos((
real(jy)*dy+ylat0)*pi180)
443 vi=vv(ix,jy,iz,n)*dyconst
446 if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
447 (height(iz).le.uvwzlev(ix,jy,kz)))
then
448 dz1=height(iz)-uvwzlev(ix,jy,kz-1)
449 dz2=uvwzlev(ix,jy,kz)-height(iz)
462 dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
463 dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
464 dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
466 dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
467 dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
468 dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
470 ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)
483 do jy=int(switchnorthg)-2,nymin1
484 ylat=ylat0+
real(jy)*dy
486 xlon=xlon0+
real(ix)*dx
488 call
cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
489 vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
499 xlon=xlon0+
real(nx/2-1)*dx
501 ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
502 vv(nx/2-1,nymin1,iz,n)**2)
503 if(vv(nx/2-1,nymin1,iz,n).lt.0.)
then
504 ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
505 vv(nx/2-1,nymin1,iz,n))-xlonr
506 elseif (vv(nx/2-1,nymin1,iz,n).gt.0.)
then
507 ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
508 vv(nx/2-1,nymin1,iz,n))-xlonr
512 if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
513 if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
519 uuaux=-ffpol*sin(xlonr+ddpol)
520 vvaux=-ffpol*cos(xlonr+ddpol)
521 call
cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
526 uupol(ix,jy,iz,n)=uupolaux
527 vvpol(ix,jy,iz,n)=vvpolaux
539 wdummy=wdummy+ww(ix,jy,iz,n)
541 wdummy=wdummy/
real(nx)
544 ww(ix,jy,iz,n)=wdummy
556 do jy=0,int(switchsouthg)+3
557 ylat=ylat0+
real(jy)*dy
559 xlon=xlon0+
real(ix)*dx
561 call
cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
562 vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
571 xlon=xlon0+
real(nx/2-1)*dx
573 ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
574 vv(nx/2-1,0,iz,n)**2)
575 if(vv(nx/2-1,0,iz,n).lt.0.)
then
576 ddpol=atan(uu(nx/2-1,0,iz,n)/ &
577 vv(nx/2-1,0,iz,n))+xlonr
578 elseif (vv(nx/2-1,0,iz,n).gt.0.)
then
579 ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
580 vv(nx/2-1,0,iz,n))-xlonr
584 if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
585 if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
591 uuaux=+ffpol*sin(xlonr-ddpol)
592 vvaux=-ffpol*cos(xlonr-ddpol)
593 call
cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
598 uupol(ix,jy,iz,n)=uupolaux
599 vvpol(ix,jy,iz,n)=vvpolaux
611 wdummy=wdummy+ww(ix,jy,iz,n)
613 wdummy=wdummy/
real(nx)
616 ww(ix,jy,iz,n)=wdummy
628 lsp=lsprec(ix,jy,1,n)
629 convp=convprec(ix,jy,1,n)
633 pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
634 rh=qv(ix,jy,kz,n)/
f_qvsat(pressure,tt(ix,jy,kz,n))
637 if ((lsp.gt.0.01).or.(convp.gt.0.01))
then
639 cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
640 height(kz)-height(kz-1)
641 if (lsp.ge.convp)
then
650 if (rain_cloud_above.eq.1)
then
651 if (lsp.ge.convp)
then
subroutine, public cc2gll(strcmp, xlat, xlong, ue, vn, ug, vg)
real function f_qvsat(p, t)