subroutine vinterp(prin,h,mr,u,v,nzin,pd,t,mreta,uetat,vetat,im,jm,nm,lm) use dynam implicit none integer, intent(in) :: im,jm,nm,lm !!! include 'dynam_comm.h' include 'const.h' integer::nzin,nzinm1,ld,l,l1,l2,l3,i,j,n,ivi,lp1 real,dimension(nzin)::prin real,dimension(0:im+1,0:jm+1,nm,nzin)::h,u,v,mr real,dimension(0:im+1,0:jm+1,nm)::pd real,dimension(0:im+1,0:jm+1,nm,lm)::t,mreta,uetat,vetat integer,dimension(lm)::ldt1 real,dimension(lm)::dlt real::prinm,rprin,h1b,h2,h3,alpgt,hetaij,alpeta & ,tetamax,alpetk,uld,vld,cf,ukb,vkb,phub,phlb,qseta, & rheta real,dimension(nzin)::alpr,alp,alpsq real,allocatable,dimension(:,:,:,:)::heta real,allocatable,dimension(:,:,:)::pdvp !!!!!!!!! ! real,dimension(lm) :: deta, aeta ! real,dimension(lp1) :: eta !!!!!!!!! integer,dimension(0:im+1,0:jm+1,nm)::lhgt real,dimension(0:im+1,0:jm+1,nm)::hgt,ref real,parameter::g=9.80616,r=287.04,pq0=379.90516,a2=17.2693882, & a3=273.16,a4=35.86,tresh=0.95 real::alp1,alp2,pres !!! lp1=lm+1 print *,'vinterp start' open(1,file='forVt.dat',status='old',form='unformatted') read(1)lhgt,hgt,ref close(1) call etaPr(prin,eta,pt,ldt1,nzin,lm) print *,'Vertically interpolate data to ETA grid.' print *,'input levels, output levels', nzin,LM print *,'in interp....pt= ', pt ! !c !----------------------------------------------------------------------- !c !c *** Vertical interpolation: quadratic interpolation of !c heights (equivalent to temperature being linear in ln p) !c and linear interpolation of winds. !c !c *** Computation of the 'sea level' pressure difference. !c nzinm1=nzin-1 prinm=prin(nzin) rprin=prinm/prin(nzinm1) do ld=1,nzinm1 alpr(ld)=alog(prin(ld+1)/prin(ld)) enddo do ld=1,nzin alp(ld)=alog(prin(ld)) alpsq(ld)=alp(ld)**2 enddo do l=1,lm l1=ldt1(l) dlt(l)=alpr(l1+1)*(alp(l1+2)-alp(l1))*(-alpr(l1)) enddo !Cmp *************** START FIRST SET OF LOOPS ******************** do n=1,nm do j=1,jm do i=1,im l=min(lm,lhgt(i,j,n)) l1=ldt1(l) l2=l1+1 l3=l1+2 h1b=h(i,j,n,l1) h2 =h(i,j,n,l2) h3 =h(i,j,n,l3) call quadInterpLnPr(h1b,h2,h3,alp(l1),alp(l2),alp(l3) & ,hgt(i,j,n),alpgt) pd(i,j,n)=(exp(alpgt)-pt)/ref(i,j,n) enddo enddo enddo ALLOCATE(HETA(0:IM+1,0:JM+1,nm,LM)) do l=lm,1,-1 do n=1,nm do i=1,im do j=1,jm l1=ldt1(l) l2=l1+1 l3=l1+2 h1b=h(i,j,n,l1) h2 =h(i,j,n,l2) h3 =h(i,j,n,l3) alpeta=alog(pt+pd(i,j,n)*eta(l)) call quadInterpH(h1b,h2,h3,alp(l1),alp(l2),alp(l3) & ,hetaij,alpeta) heta(i,j,n,l)=hetaij h1b=mr(i,j,n,l1) h2 =mr(i,j,n,l2) h3 =mr(i,j,n,l3) alpeta=alog(pt+pd(i,j,n)*aeta(l)) call quadInterpH(h1b,h2,h3,alp(l1),alp(l2),alp(l3) & ,hetaij,alpeta) mreta(i,j,n,l)=hetaij enddo enddo enddo enddo ALLOCATE(PDVP(0:IM+1,0:JM+1,nm)) do n=1,nm DO J=1,JM-1 DO I=1,IM-1 PDVP(I,J,n)=0.25*(PD(I,J,n)+PD(I+1,J,n) & +PD(I,J+1,n)+PD(I+1,J+1,n)) ENDDO ENDDO enddo !c *** Ground surface heights converted to geopotentials, and !c geopotential to temperature inversion. !c tetamax = 0 do n=1,nm do j=1,jm do i=1,im phub=0. do ivi=1,lm l=lm+1-ivi phlb=phub phub=g*heta(i,j,n,l) !JBF: ---------------------------------------------------- ! write(*,*) "l=", l, "t(i,j,n,l)=", t(i,j,n,l) !! if (t(i,j,n,l) > 320 ) write(200,*) "l=", l, "t(i,j,n,l)=", t(i,j,n,l) , "qseta=", qseta !!!CHOU t(i,j,n,l)=-(phlb-phub)*(pt+aeta(l)*pd(i,j,n)) & !!!CHOU /(r*deta(l)*pd(i,j,n)) !!?? !JBF: !! if (t(i,j,n,l) > 320 ) write(200,*) "l=", l, "i=",i, "j=",j, "t(i,j,n,l)=", t(i,j,n,l) , "qseta=", qseta !!!! if (t(i,j,n,l) > 320 ) write(200,*) "l=", l, "n=", n, "i=",i, "j=",j, "aeta(l)=", aeta(l), "deta(l)=", deta(l), "pd(i,j,n)=", pd(i,j,n), "heta(i,j,n,l)=", heta(i,j,n,l),"phlb=", phlb,"phub=",phub, "t(i,j,n,l)=", t(i,j,n,l) !!!!!!!!!!!!1 alp1=alog(pt+pd(i,j,n)*eta(l+1)) alp2=alog(pt+pd(i,j,n)*eta(l)) t(i,j,n,l)=-(phlb-phub)/(r*(alp1-alp2)) !CHOU if ( l == lm ) t(i,j,n,l) = 200/(r*(alp1-alp2)) !CHOU ! if (t(i,j,n,l) > 320 ) write(210,*) "l=", l, "n=", n, "i=",i, "j=",j, "eta(l)=", eta(l), "pd(i,j,n)=", pd(i,j,n), "heta(i,j,n,l)=", heta(i,j,n,l),"phlb=", phlb,"phub=",phub, "t(i,j,n,l)=", t(i,j,n,l), "alp1=",alp1, "alp2=", alp2 !GSM if (i == 123 .and. j == 377 .and. n == 2) write(210,*) "l=", l, "n=", n, "i=",i, "j=",j, "eta(l)=", eta(l),"aeta(l)=", aeta(l), & !GSM & "deta(l)=", deta(l), "pd(i,j,n)=", pd(i,j,n), "heta(i,j,n,l)=", heta(i,j,n,l),"phlb=", phlb,"phub=",phub, "t(i,j,n,l)=", t(i,j,n,l), "alp1=",alp1, "alp2=", alp2, "phlb-phub=", phlb-phub, "hgt(i,j,n)=", hgt(i,j,n) if (t(i,j,n,l).gt.tetamax) tetamax=t(i,j,n,l) !!!!!!!!!!!!teste ANDRE E JOAO !Cmp substitute the TETA values at lowest level? ! if (ivi .eq. lm .and. t(i,j,n,lm) > 320) then ! t(i,j,n,lm)=t(i,j,n,lm-1) ! endif ! if (t(i,j,n,l) > 320 ) write(200,*) "l=", l, "i=",i, "j=",j, "pd(i,j,n)=", pd(i,j,n), "heta(i,j,n,l)=", heta(i,j,n,l), "t(i,j,n,l)=", t(i,j,n,l) !!!!!!!!!!!teste ANDRE E JOAO t(i,j,n,l)=amin1(t(i,j,n,l),325.) t(i,j,n,l)=amax1(t(i,j,n,l),150.) ! if (t(i,j,n,l) > 320 ) write(220,*) "l=", l, "n=", n, "i=",i, "j=",j, "eta(l)=", eta(l), "pd(i,j,n)=", pd(i,j,n), "heta(i,j,n,l)=", heta(i,j,n,l),"phlb=", phlb,"phub=",phub, "t(i,j,n,l)=", t(i,j,n,l), "alp1=",alp1, "alp2=", alp2 !GSM if (i == 123 .and. j == 377 .and. n == 2) write(220,*) "l=", l, "n=", n, "i=",i, "j=",j, "eta(l)=", eta(l),"aeta(l)=", aeta(l), & !GSM & "deta(l)=", deta(l), "pd(i,j,n)=", pd(i,j,n), "heta(i,j,n,l)=", heta(i,j,n,l),"phlb=", phlb,"phub=",phub, "t(i,j,n,l)=", t(i,j,n,l), "alp1=",alp1, "alp2=", alp2, "phlb-phub=", phlb-phub, "hgt(i,j,n)=", hgt(i,j,n) !JBF ! if (t(i,j,n,1) .gt. 248.3767) then ! t(i,j,n,1)=248.3767 ! endif if (eta(l+1) .le. ref(i,j,n)) then ! pres=exp(0.5*(alp1+alp2)) qseta=pq0/(pd(i,j,n)*aeta(l)+pt) & ! qseta=pq0/(pres) & *exp(a2*(t(i,j,n,l)-a3)/(t(i,j,n,l)-a4)) !JBF: ---------------------------------------------------- ! if (t(i,j,n,l) > 320 ) write(200,*) "l=", l, "t(i,j,n,l)=", t(i,j,n,l) , "qseta=", qseta ! rheta=mreta(i,j,n,l)/qseta endif rheta=min(rheta,tresh) !!!!!!!!!! mreta(i,j,n,l)=rheta*qseta !!?? mreta(i,j,n,l)=max(0.,mreta(i,j,n,l)) !!!!!!!!!! t(i,j,n,l)=t(i,j,n,l)/(mreta(i,j,n,l)*0.61+1) !!?? enddo !c !C************************ !c ************ Now redefine pd to have it equal to ps-pt. !C************************** !c pd(i,j,n)=ref(i,j,n)*pd(i,j,n) enddo enddo enddo !JBF !!! close(200) !JBF DEALLOCATE(HETA) DO 400 L=1,LM DO n=1,nm do j=1,jm do i=1,im ALPETK=ALOG(PT+PDVP(I,J,n)*AETA(L)) ld=2 do while (ALPETK .GT. ALP(LD) .AND. LD .LT. nzin) ld=ld+1 enddo ULD=U(I,J,n,LD) VLD=V(I,J,n,LD) CF=(ALP(LD)-ALPETK)/ALPR(LD-1) UKB=ULD+(U(I,J,n,LD-1)-ULD)*CF VKB=VLD+(V(I,J,n,LD-1)-VLD)*CF UETAT(I,J,n,L)=UKB VETAT(I,J,n,L)=VKB ENDDO enddo enddo 400 CONTINUE DEALLOCATE(PDVP) !c ! *** Set u, v equal to zero at points below the ground. !c do l=1,lm do n=1,nm do j=1,jm-1 do i=1,im-1 if (eta(l+1) .gt. ref(i,j,n)) then uetat(i-1,j-1,n,l)=0. vetat(i-1,j-1,n,l)=0. uetat(i-1,j,n,l)=0. vetat(i-1,j,n,l)=0. uetat(i,j-1,n,l)=0. vetat(i,j-1,n,l)=0. uetat(i,j,n,l)=0. vetat(i,j,n,l)=0. endif enddo enddo enddo enddo print *,"finish vinterp" END SUBROUTINE vinterp