subroutine eta_interp(nnn) c implicit none include 'ecommons.h' c integer nxin,nyin,nzin . ,jday,idat(5) . ,i,l,n,ll,NN,nsfcfld,j,nnn c integer gds(200) real dummy C real,allocatable:: htin(:,:,:) . ,mrin(:,:,:),uwin(:,:,:),vwin(:,:,:),sfcgrid(:,:,:) . ,tpin(:,:,:) Chou Temperature included real prin(nzin),dtr c real hthi(im,jm,nzin) . ,mrhi(im,jm,nzin) . ,uwhi(im,jm,nzin) . ,vwhi(im,jm,nzin) . ,sfcgridhi(im,jm,12),qmn,qmx . ,tphi(im,jm,nzin) Chou real,allocatable:: hgt(:,:) real sm(im,jm) Chou integer k parameter(dtr=.017453293) c character*2 gproj character*100 dummyChar c logical first data first/.true./ c common /sectorsize/nxin,nyin,nzin c save idat c_______________________________________________________________________________ c c write(6,*) nzin,im,jm c c *** Get time increment for boundary conditions. c if (first) then idat(5)=-nint(tboco) first=.false. c print *,'first:',idat(5) endif c idat(1)=imonth idat(2)=idate idat(3)=iyear idat(4)=istrtim idat(5)=idat(5)+nint(tboco) c c *** Read in input upper air gridded data set. c call get_sector_size(init_in(nnn),init_gdsdir) print*,'Data Dimension',NXIN,NYIN,NZIN ALLOCATE(HTIN(NXIN,NYIN,NZIN)) ALLOCATE(UWIN(NXIN,NYIN,NZIN)) ALLOCATE(VWIN(NXIN,NYIN,NZIN)) ALLOCATE(MRIN(NXIN,NYIN,NZIN)) ALLOCATE(TPIN(NXIN,NYIN,NZIN)) !chou: TP ALLOCATE(SFCGRID(NXIN,NYIN,12)) call getdata(init_in(nnn),prin,uwin,vwin,htin,mrin,tpin,sfcgrid) C write(6,*) 'Q after getdata' do L=1,NZIN qmn=99999. qmx=-99999. DO J=1,NYIN do I=1,NXIN if (mrin(i,j,l) .gt. qmx) qmx=mrin(i,j,l) if (mrin(i,j,l) .lt. qmn) qmn=mrin(i,j,l) enddo enddo C write(6,*) 'raw mrin extremes: ', L,qmn,qmx enddo Chou print if (.false.) then open(unit=3,file='avnin.dat',form='unformatted', & access='sequential',recl=nxin*nyin*4) do k=nzin,1,-1 write(3) ((htin(i,j,k),i=1,nxin),j=1,nyin) enddo do k=nzin,1,-1 write(3) ((uwin(i,j,k),i=1,nxin),j=1,nyin) enddo do k=nzin,1,-1 write(3) ((vwin(i,j,k),i=1,nxin),j=1,nyin) enddo do k=nzin,1,-1 write(3) ((mrin(i,j,k),i=1,nxin),j=1,nyin) enddo do k=nzin,1,-1 write(3) ((tpin(i,j,k),i=1,nxin),j=1,nyin) enddo endif c c *** Horizontally interpolate input data from native grid to ETA grid. c init_gdsdir=init_in(nnn) n=index(init_gdsdir,' ')-1 init_gdsdir=init_gdsdir(1:n)//'.gdsinfo' write(6,*) 'trying to open ', init_gdsdir(1:n) n=index(init_gdsdir,' ')-1 open(14,file=init_gdsdir(1:n),form='formatted', + status='old') rewind 14 read(14,*) NZIN read(14,*) dummyChar do L=1,NZIN read(14,*) dummy enddo gds(1:200)=0 read(14,*) gds(1) read(14,*) gds(2) read(14,*) gds(3) read(14,*) gds(4) read(14,*) gds(5) read(14,*) gds(6) read(14,*) gds(7) read(14,*) gds(8) read(14,*) gds(9) read(14,*) gds(10) read(14,*) gds(11) read(14,*) gds(20) write(6,*) 'gds(1-14): ', (gds(I),I=1,14) call hinterp(htin,gproj,hthi,nzin,1,gds) DEALLOCATE(HTIN) call hinterp(mrin,gproj,mrhi,nzin,2,gds) DEALLOCATE(MRIN) call hinterp(tpin,gproj,tphi,nzin,1,gds) !Chou inc TP DEALLOCATE(TPIN) !Chou inc TP write(6,*) 'calling wind interp' call hinterp_wind(uwin,vwin,gproj,uwhi,vwhi,nzin,5,gds) DEALLOCATE(UWIN) DEALLOCATE(VWIN) call hinterp(sfcgrid,gproj,sfcgridhi,12,6,gds) C ALLOCATE(HGT(IM,JM)) open(1,file=topo_out,status='old',form='unformatted') rewind(1) read(1) hgt,sm close(1) DEALLOCATE(HGT) C if (.not. REAN .and. GRIBSOIL .and. idat(5) .eq. 0) then write(6,*) 'calling coastal!!!!!!!' call coastalfix(sfcgridhi,sfcgrid,sm,im,jm,nxin,nyin,dphd,dlmd, + wbd,sbd,tph0d,tlm0d,gds) write(6,*) 'return coastal' endif DEALLOCATE(SFCGRID) c c *** Vertically interpolate input data from native grid to ETA grid. c *** And write output files. c Chou inc TP call vinterp(prin,hthi,mrhi,uwhi,vwhi,tphi,sfcgridhi,nzin,idat) c end subroutine eta_interp C ********************* END MAIN ROUTINE ************************************ C C subroutine hinterp(fieldin,gproj,fieldint,nzmax,mark,gds) c implicit none c include 'ecommons.h' c c *** Input data on native grid. c integer nxsec,nysec,nzsec,mark,L,nzmax,n c real fieldin(nxsec,nysec,nzsec) c c *** Input data on horizontal eta grid, native vertical grid. c real fieldint(im,jm,nzsec) REAL,ALLOCATABLE:: tmpout(:,:),tmp(:,:) c real*4 coh(3,im,jm),inh(4,im,jm),jnh(4,im,jm) . ,cov(3,im,jm),inv(4,im,jm),jnv(4,im,jm) . ,ald(im,jm),apd(im,jm) integer i,j,k integer GDS(200) character*2 gproj c common /sectorsize/nxsec,nysec,nzsec c print *,'Horizontally interpolate data to ETA grid.... ', mark c if (gds(1) .eq. 3) gproj='LC' if (gds(1) .eq. 5) gproj='PS' if (gds(1) .eq. 0) gproj='LL' write(6,*) 'calling with gproj= ', gproj call gtll(coh,inh,jnh,cov,inv,jnv,ald,apd,nxsec,nysec,gproj,gds) c c *** Horizontally interpolate input data (ht, mr, u, and v) to ETA grid c using bilinear interpolation. c DO L=1,nzmax ALLOCATE (TMP(NXSEC,NYSEC)) tmp(1:nxsec,1:nysec)=fieldin(1:nxsec,1:nysec,l) ALLOCATE (TMPOUT(IM,JM)) if (mark .ne. 4 .and. mark .ne. 5) then call bilinb(coh,inh,jnh,im,jm,nxsec,nysec,tmp,tmpout) else C write(6,*) 'interpolating winds ' call bilinb(cov,inv,jnv,im,jm,nxsec,nysec,tmp,tmpout) endif DEALLOCATE(TMP) fieldint(1:im,1:jm,l)=tmpout(1:im,1:jm) DEALLOCATE(TMPOUT) ENDDO c c *** Check that mixing ratio is not less than zero. c if (mark .eq. 2) then write(6,*) 'checking mixing ratio for neg values' where(fieldint<=0.) fieldint=0. endif c return end c C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C subroutine hinterp_wind(uin,vin,gproj,uint,vint,nzmax,mark,gds) c implicit none c include 'ecommons.h' c c *** Input data on native grid. c integer nxsec,nysec,nzsec,mark,L,nzmax,n c real uin(nxsec,nysec,nzsec),vin(nxsec,nysec,nzsec) c c *** Input data on horizontal eta grid, native vertical grid. c real uint(im,jm,nzsec),vint(im,jm,nzsec) REAL,ALLOCATABLE:: utmpout(:,:),utmp(:,:),utmpout1(:,:) REAL,ALLOCATABLE:: vtmpout(:,:),vtmp(:,:),vtmpout1(:,:) c real*4 coh(3,im,jm),inh(4,im,jm),jnh(4,im,jm) . ,cov(3,im,jm),inv(4,im,jm),jnv(4,im,jm) . ,ald(im,jm),apd(im,jm),dtr parameter(dtr=.017453293) integer i,j,k integer GDS(200) character*2 gproj c common /sectorsize/nxsec,nysec,nzsec c print *,'Horizontally interpolate wind data to ETA grid.... ' c if (gds(1) .eq. 3) gproj='LC' if (gds(1) .eq. 5) gproj='PS' if (gds(1) .eq. 0) gproj='LL' write(6,*) 'calling with gproj= ', gproj call gtll(coh,inh,jnh,cov,inv,jnv,ald,apd,nxsec,nysec,gproj,gds) c c *** Horizontally interpolate input data (ht, mr, u, and v) to ETA grid c using bilinear interpolation. c DO L=1,nzmax ALLOCATE (UTMP(NXSEC,NYSEC)) ALLOCATE (VTMP(NXSEC,NYSEC)) Cold DO N=1,2 do j=1,nysec do I=1,nxsec utmp(i,j)=uin(i,j,l) vtmp(i,j)=vin(i,j,l) enddo enddo ALLOCATE (UTMPOUT(IM,JM)) ALLOCATE (VTMPOUT(IM,JM)) ALLOCATE (UTMPOUT1(IM,JM)) ALLOCATE (VTMPOUT1(IM,JM)) call bilinb(cov,inv,jnv,im,jm,nxsec,nysec,utmp,utmpout) call bilinb(cov,inv,jnv,im,jm,nxsec,nysec,vtmp,vtmpout) DEALLOCATE(UTMP) DEALLOCATE(VTMP) call ltlwin(ald,apd,utmpout,vtmpout,dtr,utmpout1,vtmpout1) do j=1,jm do I=1,im uint(i,j,l)=utmpout1(i,j) vint(i,j,l)=vtmpout1(i,j) enddo enddo Cold ENDDO DEALLOCATE(UTMPOUT) DEALLOCATE(VTMPOUT) DEALLOCATE(UTMPOUT1) DEALLOCATE(VTMPOUT1) ENDDO c return end c========================================================================== C subroutine vinterp(spl,h,q,u,v,t,sfcgrid,ldm,idat) c c *** Code to vertically interpolate gridded init data from c native vertical grid to ETA vertical grid. c c *** Original code received from the U. of Athens and modfied at FSL. c c *** Setting up of the vertical grid variables, and c vertical interpolation, pressure to eta surfaces c Fedor Mesinger c implicit none c include 'ecommons.h' include 'econstants.h' c real*4 a2,a3,a4,pq0 real*4 a12,ai2,ai4,vtmpc1,vtmpc2,corr parameter (a2=17.2693882E0,a3=273.16E0,a4=35.86E0,pq0=379.90516) !Chou Consider below zero temperature parameter (ai2=21.875,ai4=7.66,vtmpc1=.60776868,vtmpc2=.85930818) !Chou C integer*4 ldm,negsum,ll,NN c C ldm REAL H(IM,JM,LDM),U(IM,JM,LDM),V(IM,JM,LDM),Q(IM,JM,LDM) REAL T(IM,JM,LDM) !chou inc Tp C lm REAL, ALLOCATABLE:: HETA(:,:,:),QETA(:,:,:),UETA(:,:,:) REAL, ALLOCATABLE:: VETA(:,:,:),TETA(:,:,:),PDVP(:,:),PHIS(:,:) REAL, ALLOCATABLE:: HGT(:,:),SM(:,:),PD(:,:) real*4 . ext(lm,2) . ,hetaij . ,sfcgrid(im,jm,12) . ,pdrv,petatmp real uetaij(lm),vetaij(lm) . ,href(ldm),eta(lmp1),peta(lmp1) . ,pdb(kb,2),tb(kb,lm,2),ub(kb,lm,2),vb(kb,lm,2) . ,qb(kb,lm,2) . ,etal(lm),zeta(lmp1),ref(im,jm) . ,hfp(2,2),hfpa(im,jm),hrsd(im,jm),dz(im,jm) . ,etas(im,jm),wm(im,jm) . ,seta(4),deta(lm),dum(lm) . ,alpq,bq . ,q3,q3m1,q3m2,hldm,cq,rheta,tresh . ,tetam,cf,uld,vld,qseta,tetamax . ,phub,phlb,q1,q2b,prfs,tmj2,gorg,rgog . ,pt,dpd,dpu,x,h1b,h2,h3,rspl,splm,gfac,alpt . ,t1,t2,t3,t3m1,t3m2,bt,ct,zdif,adjfac,tsavetmp,t2b, . ukb,vkb . ,tkb,qkb, ttld, qqld, rseta !Chou 20240929 vinterp T Q as U V Cmp real(8) alpeta,alpetai(im,jm), alp(ldm), alpr(ldm+1),alpsq(ldm) real(8) dlt(lm),alpub,alplb,b,c,alpgt,h3m1,h3m2,alpetk c integer*4 lhgt(im,jm),lnew(im,jm),ldt1(lm),ldum(lm),MaskInt(lm) . ,idat(5),lmin(im,jm) . ,mi2,ip1,jp1,jm1b,ldm1,ivi . ,is,js,if,jf,l1,l2,l3,ldmm1,nk . ,kt,lin,ii,iii . ,lnbmx,lmnn,lmne,lmns,lmnw,ld,lij,lrais . ,mmx,ntsd,i,j,l,m,n,ihw(jm),ihe(jm),ivw(jm),ive(jm) . ,ivl,ivh,ihl,ihh,mnknt,idxx(im*jm),jdxx(im*jm),numk . ,kntx,idx(im*jm),jdx(im*jm),kntin,kp1 real wt(LM),wb(LM), DPSLP integer it_lbc(LM), Llbc c character*8 fname c logical wndls,hmask(im,jm),vmask(im,jm),run c real*4 pdmin,pdmax, spl(ldm) c_______________________________________________________________________________ c print *,'Vertically interpolate data to ETA grid.' write(6,*) 'input levels, output levels', LDM,LM negsum=0 ldmm1=ldm-1 !Chou 20240929 tresh=0.97 ! Chou increased threshold tresh=0.999 ! Chou 20240929 increased threshold from 0.97 to 0.999 c kt=idat(5) Cmp pt=ptinp Cmp write(6,*) 'in interp....pt= ', pt write(6,*) 'in interp....ptinp= ', ptinp Czj !-------------indirect indices for tom's set-up------------------------- do j=1,jm ihw(j)=-1+mod(j+1,2) ihe(j)=ihw(j)+1 ivw(j)=-1+mod(j,2) ive(j)=ivw(j)+1 enddo !---------------------------------------------------------------------- Czj c c----------------------------------------------------------------------- c print=.true. c mountains to be four-point ('regular' or 'silhouette') c averaged, so as to have the same ground topography at groups c of four neighboring height points, or not four-point averaged fpmnts=.false. c siluet=.true. c hour00 fields to be created for a sigma mode/not a sigma mode c experiment c sigma=.false. c if a sigma mode experiment, mountains to be removed 6 grids c lines along the boundary/not to be removed c mrmsxl=.false. cda mrmsxl=.not.sigma c no mountains experiment/not a no mountains experiment c nomnts=.false. c c *** Define some constants needed for vertical interpolation. c gorg=g/(r*gamma) rgog=r*gamma/g alpt=log(pt) !Chou changed to log c------- eta at the interfaces of the model layers ------------------ C read in Eta levels from file open(unit=16,file='deta',form='unformatted' + ,access='sequential') REWIND 16 READ(16)DETA,LDUM eta(1)=0.0 DO L=1,LM ETA(L+1) = ETA(L)+DETA(L) c *** Compute half eta values. ETAL(L) = 0.5*(ETA(L)+ETA(L+1)) peta(l)=(prf0-pt)*eta(l)+pt zeta(l)=T0*(1.-((pt+eta(l) & *(prf0-pt))/prf0)**rgog)/gamma Cnew C write(6,*) 'peta(l), spl(i), i ', peta(l),spl(i),i ENDDO ETA(LMP1)=1.0 CCHOU 20260901 C Construct weights for vertical interpolation. If Peta - SPL <= 5hpa C then Eta_var = GCM_var : no interpolation Llbc=1 DOLM:Do L=1,LM do while (peta(l) .ge. spl(Llbc) .and. Llbc .lt. ldm) !CHOU LDM=23 EC Llbc=Llbc+1 enddo IF (ABS(PETA(L)-SPL(Llbc-1)) .le. 500.) then wt(l)=1. wb(l)=0. it_lbc(l)=Llbc continue !! enddo dolm ELSEIF (ABS(PETA(L)-SPL(Llbc)).le. 500.) then wt(l)=0. wb(l)=1. it_lbc(l)=Llbc continue !! enddo dolm ELSE DPSLP=(SPL(Llbc)-SPL(Llbc-1)) wt(l)=(SPL(Llbc)-PETA(l))/DPSLP wb(l)=(PETA(l)-SPL(Llbc-1))/DPSLP it_lbc(l)=Llbc continue !! enddo dolm !GSM ENDIF ENDIF ENDDO DOLM c CChou retirou esta parte e substituiu com o trecho de cima c *** Compute ldt1. c C do l=1,lm C i=1 C ii=1 C do while (peta(l) .ge. spl(i) .and. i .lt. ldm) !CHOU LDM=23 EC C i=i+1 C enddo C if (.not.lInterpAllLev) then C do ii=1, ldm C if ((abs(peta(l)-spl(ii)))*1E-2 .lt.2 )then C MaskInt(l)=1+ii C iii=ii C exit C endif C MaskInt(l)=0 C iii=i-1 C enddo C else C MaskInt(l)=0 C endif CC C if (i .le. 2) then C ldt1(l)=1 C elseif (i .eq. ldm) then C ldt1(l)=ldm-2 C else C dpu=abs(spl(i-2)-peta(l)) C dpd=abs(spl(i+1)-peta(l)) C if (dpu .gt. dpd) then C ldt1(l)=i-1 C else C ldt1(l)=i-2 C endif C endif C enddo dolm c c c c *** Compute eta level heights. c write(6,*) 'lowest zeta values' do l=1,lmp1 zeta(l)=t0*(1.-((pt+eta(l)*(prf0-pt))/prf0)**rgog)/gamma C if (l .ge. lmp1-10) write(6,*) 'zeta(l)... ', l, zeta(l) enddo c c *** Read topo heights and sea masks. C C This version expects hgt, sm to be (im,jm) and NOT (imt,jmt) C make sure etatopo.f handles this properly!!!! C ALLOCATE(HGT(IM,JM)) ALLOCATE(SM(IM,JM)) open(1,file=topo_out,status='old',form='unformatted') rewind(1) read(1) hgt,sm close(1) c *** Zero out heights, if no mountains option has been selected. c C if (nomnts) then C do j=1,jm C do i=1,im C if (hgt(i,j).gt.0.) hgt(i,j)=0.004 C ref(i,j)=1. C enddo C enddo C goto 922 C endif c doout902: do j=1,jm doin902: do i=1,im prfs=prf0*((t0-gamma*hgt(i,j))/t0)**gorg etas(i,j)=(prfs-pt)/(prf0-pt) C C Cmp through this point hgt represents the unchanged raw hgt value C do903: do ivi=1,lm l=lmp1-ivi if (etas(i,j).lt.eta(l)) cycle do903 cvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv if (etas(i,j).lt.etal(l)) then c if (etas(i,j).lt.0.10*(9.*etal(l)+eta(l-1))) then caaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa ref (i,j)= eta(l) hgt (i,j)=zeta(l) lhgt(i,j)= l else ref (i,j)=eta(l+1) if (l+1 .eq. lmp1 .and. hgt(i,j) .gt. 0.) then hgt(i,j)=0.004 else hgt(i,j)=zeta(l+1) endif lhgt(i,j)= l+1 endif c cycle doin902 end do do903 end do doin902 end do doout902 c c c if(.not.fpmnts): c find points having no neighboring wind point because of mountains c on all four sides (sunken points), and c * if land, just raise to reach the lowest wind; c * if an isolated water point, or a water point with only one c neighboring water point and at sea level: raise, and declare c land (in order not to trigger convection over artificially c elevated warm water); c * otherwise: c remove (level off) land at one of the four neighboring wind c points -- the one having the smallest three-point averaged c elevation -- so as to restore one of the wind points c c n.b.: i am in the "old" i,j indexing, with i being incremented c in each column c if (.not. fpmnts) then c c Save first the sea mask, sm, into a back-up mask ("water mask", c wm) to use in testing the neighbors, since sm will be changed c along the way; and initialize lmin, level beyond which points c should not be raised because they have participated in leveling c off at that level, and lnew c do j=1,jm do i=1,im wm (i,j)=sm (i,j) lmin(i,j)=2 lnew(i,j)=lhgt(i,j) enddo enddo c do j=4,jm-3 ihl=2+mod(j,2) ihh=im-2 do i=ihl,ihh c c ************ Find the maximum elevation (minimum l) of four neighboring wind c points -- to the west, east, south and north c lmnw=min(lhgt(i+ihw(j),j-1),lhgt(i-1,j),lhgt(i+ihw(j),j+1)) lmne=min(lhgt(i+ihe(j),j-1),lhgt(i+1,j),lhgt(i+ihe(j),j+1)) lmns=min(lhgt(i,j-2),lhgt(i+ihw(j),j-1),lhgt(i+ihe(j),j-1)) lmnn=min(lhgt(i+ihw(j),j+1),lhgt(i+ihe(j),j+1),lhgt(i,j+2)) c c ************ Which one is the lowest of the four maximum elevations? c lnbmx=max(lmnw,lmne,lmns,lmnn) c c ************ Is the points original elevation below that of the lowest c available original wind, and, at the same time, in case it c participated in leveling off is it still below the level of the c lowest wind created by leveling off? c if (lhgt(i,j) .gt. lnbmx .and.ref(i,j) & .gt. eta(lmin(i,j))) then lrais=max0(lnbmx,lmin(i,j)) if (sm(i,j) .eq. 0.) then ref (i,j)= eta(lrais) hgt (i,j)=zeta(lrais) lnew(i,j)= lrais c else if (wm(i+ihw(j),j-1)+wm(i+ihe(j),j-1) & +wm(i+ihw(j),j+1)+wm(i+ihe(j),j+1).eq.0..or. & wm(i+ihw(j),j-1)+wm(i+ihe(j),j-1) & +wm(i+ihw(j),j+1)+wm(i+ihe(j),j+1).eq.1..and. & lhgt(i,j).eq.lm+1) then ref (i,j)= eta(lrais) hgt (i,j)=zeta(lrais) lnew(i,j)= lrais c else c c *************** There is a neighboring water point, and if it is only a c single neighboring water point the water points are above sea c level. denote wind points to the east, north, west c and south by 1, 2, 3, and 4, respectively, and check the c three-point volume at each of them. c seta(1)=etas(i+ihe(j),j-1) & +etas(i+1,j)+etas(i+ihe(j),j+1) seta(2)=etas(i+ihw(j),j+1) & +etas(i+ihe(j),j+1)+etas(i,j+2) seta(3)=etas(i+ihw(j),j-1) & +etas(i-1,j)+etas(i+ihw(j),j+1) seta(4)=etas(i,j-2) & +etas(i+ihw(j),j-1)+etas(i+ihe(j),j-1) c mmx=1 do m=2,4 if(seta(m) .gt. seta(mmx)) mmx=m enddo c c ****************** n.b.: When removing now the obstructing land around the c wind point vm, lhgt is not changed. Thus, hopefully (?) c resulting topography does not depend on the direction of c the sweep within the current loop. c lij=lhgt(i,j) if (mmx .eq. 1 .or. mmx .eq. 4) then if (ref(i+ihe(j),j-1).lt.eta(lij)) then ref(i+ihe(j),j-1)=eta(lij) hgt(i+ihe(j),j-1)=zeta(lij) lnew(i+ihe(j),j-1)=lij endif lmin(i+ihe(j),j-1)=max(lmin(i+ihe(j),j-1),lij) endif c if (mmx .eq. 1) then if (ref(i+1,j) .lt. eta(lij)) then ref(i+1,j )= eta(lij) hgt(i+1,j )= zeta(lij) lnew(i+1,j )= lij endif lmin(i+1,j)=max(lmin(i+1,j),lij) endif c if (mmx .eq. 1 .or. mmx .eq. 2) then if (ref(i+ihe(j),j+1).lt.eta(lij)) then ref(i+ihe(j),j+1)=eta(lij) hgt(i+ihe(j),j+1)=zeta(lij) lnew(i+ihe(j),j+1)=lij endif lmin(i+ihe(j),j+1)=max(lmin(i+ihe(j),j+1),lij) endif c if (ref(i,j+1) .lt. eta(lij)) then ref(i ,j+2 )= eta(lij) hgt(i ,j+2 )= zeta(lij) lnew(i ,j+2 )= lij endif lmin(i,j+2)=max(lmin(i,j+2),lij) endif c if (mmx .eq. 2 .or. mmx .eq. 3) then if (ref(i+ihw(j),j+1).lt.eta(lij)) then ref(i+ihw(j),j+1)=eta(lij) hgt(i+ihw(j),j+1)=zeta(lij) lnew(i+ihw(j),j+1)=lij endif lmin(i+ihw(j),j+1)=max(lmin(i+ihw(j),j+1),lij) endif c if (mmx .eq. 3) then if (ref(i-1,j) .lt. eta(lij)) then ref(i-1,j)= eta(lij) hgt(i-1,j)= zeta(lij) lnew(i-1,j)= lij endif lmin(i-1,j)=max(lmin(i-1,j),lij) endif c if (mmx .eq. 3 .or. mmx .eq. 4) then if (ref(i+ihw(j),j-1).lt.eta(lij)) then ref(i+ihw(j),j-1)=eta(lij) hgt(i+ihw(j),j-1)=zeta(lij) lnew(i+ihw(j),j-1)=lij endif lmin(i+ihw(j),j-1)=max(lmin(i+ihw(j),j-1),lij) endif c if (mmx .eq. 4) then if (ref(i,j-2).lt.eta(lij)) then ref(i,j-2)=eta(lij) hgt(i,j-2)=zeta(lij) lnew(i,j-2)=lij endif lmin(i,j-2)=max(lmin(i,j-2),lij) endif endif enddo enddo c c ****** Any windless points left or inadvertently created? c wndls=.false. do j=4,jm-3 ihl=2+mod(j,2) ihh=im-2 ! do i=ihl,ihh ! c ************ Find the maximum elevation (minimum l) of four neighboring wind c points -- to the west, east, south and north c lmnw=min(lnew(i+ihw(j),j-1),lnew(i-1,j),lnew(i+ihw(j),j+1)) lmne=min(lnew(i+ihe(j),j-1),lnew(i+1,j),lnew(i+ihe(j),j+1)) lmns=min(lnew(i,j-2),lnew(i+ihw(j),j-1),lnew(i+ihe(j),j-1)) lmnn=min(lnew(i+ihw(j),j+1),lnew(i+ihe(j),j+1),lnew(i,j+2)) c c ************ Which one is the lowest of the four maximum elevations? c lnbmx=max(lmnw,lmne,lmns,lmnn) c c ************ Is the points new elevation below that of the lowest available c new wind? c if (lnew(i,j) .gt. lnbmx) wndls=.true. enddo enddo c endif c if (sigma) then do j=1,jm do i=1,im ref(i,j)=1. enddo enddo endif c c *** Masks. do j=1,jm do i=1,im hmask(i,j)=.false. vmask(i,j)=.false. enddo enddo !----------------------------------------------------------------------- do i=1,im hmask(i,1)=.true. hmask(i,jm)=.true. enddo do j=1,jm,2 hmask(1,j)=.true. hmask(im,j)=.true. enddo !----------------------------------------------------------------------- do i=1,im vmask(i,1)=.true. vmask(i,jm)=.true. enddo do j=2,jm-1,2 vmask(1,j)=.true. vmask(im,j)=.true. enddo !----------------------------------------------------------------------- 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 splm=spl(ldm) rspl=splm/spl(ldmm1) c Cmp *************** START FIRST SET OF LOOPS ******************** ALLOCATE(PD(IM,JM)) doout1: do j=1,jm doin1: do i=1,im Cmp approach started 2/5/99....use lhgt to only use the closest Cmp to ground surface data...not subground data Chou alpgt=(-b-sqrt(b**2-4.*c Chou . *(-c*alpsq(l3)-b*alp(l3)+h3-hgt(i,j))))/(2.*c) C pd for "normal" points Cmp want to save the pd value to an array so have when do rest Cmp of calculations below Chou pd(i,j)=(exp(alpgt)-pt)/ref(i,j) Chou goto 533 532 hldm=h(i,j,ldm) Cmp vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv pd(i,j)=(splm*rspl**((hldm-hgt(i,j)) . /(h(i,j,ldmm1)-hldm))-pt) . /ref(i,j) Cmp ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Chou write(6,*) 'special case.... ', i,j,pd(i,j) Chou write(6,969) hldm,hgt(i,j),h(i,j,ldmm1),ref(i,j) 969 format(' hldm,hgt,h,ref ',3(f8.2,x),f7.5) write(6,*) '.........................................' Cmp *** when pd(i,j) is large, alpeta is larger than normal Cmp *** then when take differences with alp() and alpsq() which are functions Cmp *** of the input data set pressure levels get values of hetaij which do Cmp *** not match well with surrounding regions 533 continue enddo doin1 enddo doout1 C************************************************************* C NEW STUFF **************************** C************************************************************* ALLOCATE(HETA(IM,JM,LM)) ALLOCATE(QETA(IM,JM,LM)) ALLOCATE(TETA(IM,JM,LM)) !Chou in Tp 611 format('alpeta, alp(l2), b, c ',4(f14.6,1x)) C **************************************************************** C****************** end of new stuff ***************************** C***************************************************************** pdmax=0. pdmin=99999. do J=1,JM do I=1,IM if (pd(I,J)*ref(I,J).lt.pdmin .and. pd(I,J).gt.1000.) + pdmin=pd(I,J)*ref(I,J) if (pd(I,J)*ref(I,j).gt.pdmax) then pdmax=pd(I,J)*ref(I,J) endif enddo enddo write(6,*) 'in interp, pdmin (surface based)= ', pdmin write(6,*) 'in interp, pdmax (surface based)= ', pdmax c c *** Computation of winds within eta layers. c ALLOCATE(PDVP(IM,JM)) DO J=2,JM-1 DO I=2,IM-1 PDVP(I,J)=0.25*(PD(I+IVE(J),J)+PD(I+IVW(J),J) 1 +PD(I,J+1)+PD(I,J-1)) ENDDO ENDDO C*** C*** COMPUTE NORTH AND SOUTH EDGES AND RECOMPUTE EAST AND WEST SIDES C*** DO I=1,IM-1 PDVP(I,1) =0.5*(PD(I+IVE(1),1)+PD(I+IVW(1),1)) PDVP(I,JM)=0.5*(PD(I+IVE(JM),JM)+PD(I+IVW(JM),JM)) ENDDO DO J=2,JM-1,2 PDVP(1,J) =0.5*(PD(1,J+1)+PD(1,J-1)) PDVP(IM,J)=0.5*(PD(IM,J+1)+PD(IM,J-1)) ENDDO DO J=1,JM,2 PDVP(IM,J)=PDVP(IM-1,J) ENDDO DO J=3,JM-2,2 PDVP(1,J)=0.25*(PD(1+IVE(J),J)+PD(1+IVW(J),J) 1 +PD(1,J+1)+PD(1,J-1)) ENDDO C*** C*** COMPUTE WINDS (SINCE FIX IS HARDWIRED TO 0. FOR ALL K, C*** SKIP CHECKING THE VALUE OF FIX AND PROCEED) C*** ALLOCATE(UETA(IM,JM,LM)) ALLOCATE(VETA(IM,JM,LM)) write(6,*) 'past wind allocates' do400: DO L=1,LM cdir$ ivdep DO J=1,JM DO I=1,IM HETA(I,J,L) = H(I,J,it_lbc(L)-1)*wt(L) + H(I,J,it_lbc(L))*wb(L) TETA(I,J,L) = T(I,J,it_lbc(L)-1)*wt(L) + T(I,J,it_lbc(L))*wb(L) QETA(I,J,L) = Q(I,J,it_lbc(L)-1)*wt(L) + Q(I,J,it_lbc(L))*wb(L) UETA(I,J,L) = U(I,J,it_lbc(L)-1)*wt(L) + U(I,J,it_lbc(L))*wb(L) VETA(I,J,L) = V(I,J,it_lbc(L)-1)*wt(L) + V(I,J,it_lbc(L))*wb(L) teta(i,j,l)=amin1(teta(i,j,l),325.) teta(i,j,l)=amax1(teta(i,j,l),150.) qeta(i,j,l)=max(0.,qeta(i,j,l)) !JLG Esse teste melhorou o a chuva na região amazonica e ITCZ, !! porém, piorou os outros campos. Alteração descartada if (lDryT) teta(i,j,l)=teta(i,j,l)/(qeta(i,j,l)*0.61+1) !now dry T enddo enddo enddo do400 DEALLOCATE(PDVP) write(6,*) 'past PDVP deallocate' c *** Ground surface heights converted to geopotentials, and c geopotential to temperature inversion. c ALLOCATE(PHIS(IM,JM)) do j=1,jm do i=1,im phis(i,j)=g*hgt(i,j) phub=0. c c ************ Now redefine pd to have it equal to ps-pt. C************************** pd(i,j)=ref(i,j)*pd(i,j) C C**************************** C************************** enddo enddo DEALLOCATE(HETA) DEALLOCATE(HGT) c c *** Set u, v and t equal to zero at points below the ground. c Chou do j=2,jm-1 Chou ihl=1+mod(j,2) Chou ihh=im-1 ! Cmp i indices changed from zj code Chou do i=ihl,ihh Chou do l=1,lm Cmp---change this part???????????????????????????? Chou if (eta(l+1) .gt. ref(i,j)) then Chou ueta(i+ihw(j),j,l)=0. Chou veta(i+ihw(j),j,l)=0. Chou ueta(i+ihe(j),j,l)=0. Chou veta(i+ihe(J),j,l)=0. Chou ueta(i,j-1,l)=0. Chou veta(i,j-1,l)=0. Chou ueta(i,j+1,l)=0. Chou veta(i,j+1,l)=0. Chou endif Chou Chou enddo Chou enddo Chou enddo c c *** The separation of the boundary values. c n = 1 do i=1,im pdb(n,1)=pd(i, 1) pdb(n,2)=0. do l=1,lm tb(n,l,1)=teta(i, 1,l) tb(n,l,2)=0. qb(n,l,1)=qeta(i, 1,l) qb(n,l,2)=0. enddo n=n+1 enddo do i=1,im pdb(n,1)=pd(i,jm) pdb(n,2)=0. do l=1,lm tb(n,l,1)=teta(i,jm,l) tb(n,l,2)=0. qb(n,l,1)=qeta(i,jm,l) qb(n,l,2)=0. enddo n=n+1 enddo do j=3,jm-2,2 pdb(n,1)=pd( 1,j) pdb(n,2)=0. do l=1,lm tb(n,l,1)=teta( 1,j,l) tb(n,l,2)=0. qb(n,l,1)=qeta( 1,j,l) qb(n,l,2)=0. enddo n=n+1 enddo do j=3,jm-2,2 pdb(n,1)=pd(im,j) pdb(n,2)=0. do l=1,lm tb(n,l,1)=teta(im,j,l) tb(n,l,2)=0. qb(n,l,1)=qeta(im,j,l) qb(n,l,2)=0. enddo n=n+1 enddo c C SOUTH BOUNDARY - V POINTS c n=1 do i=1,im-1 do l=1,lm ub(n,l,1)=ueta(i, 1,l) ub(n,l,2)=0. vb(n,l,1)=veta(i, 1,l) vb(n,l,2)=0. enddo n=n+1 enddo C C NORTH BOUNDARY - V POINTS C do i=1,im-1 do l=1,lm ub(n,l,1)=ueta(i,jm,l) ub(n,l,2)=0. vb(n,l,1)=veta(i,jm,l) vb(n,l,2)=0. enddo n=n+1 enddo C C WEST BOUNDARY - V POINTS C do j=2,jm-1,2 do l=1,lm ub(n,l,1)=ueta( 1,j,l) ub(n,l,2)=0. vb(n,l,1)=veta( 1,j,l) vb(n,l,2)=0. enddo n=n+1 enddo C C EAST BOUNDARY - V POINTS C do j=2,jm-1,2 do l=1,lm ub(n,l,1)=ueta(im,j,l) ub(n,l,2)=0. vb(n,l,1)=veta(im,j,l) vb(n,l,2)=0. enddo n=n+1 enddo c if (kt .eq. 0) then c write(6,*) 'THE LARGEST TETA VALUE IS ', TETAMAX Cmp l=index(init_out//' ',' ')-1 if (init_out(l:l) .ne. '/') then l=l+1 init_out(l:l)='/' endif c open(1,file=init_out(1:l)//'preproc.init' . ,status='unknown',form='unformatted') rewind(1) run=.true. ntsd=0 print *,'Writing to file: ',init_out(1:l)//'preproc.init' 666 format(20(f5.0,1x)) 667 format(20(L2,1x)) write(1) run,idat(1),idat(2),idat(3),idat(4),ntsd,ueta,veta write(1) teta,qeta,pd,phis,sm,ref,eta,pt,deta,etal,zeta print *,ueta(1,1,20),ueta(11,11,20),ueta(21,21,20) print *,veta(1,1,20),veta(11,11,20),veta(21,21,20) Cmp ***** write out the surface fields write(1) sfcgrid close(1) c endif DEALLOCATE(PD) DEALLOCATE(PHIS) DEALLOCATE(SM) DEALLOCATE(TETA) DEALLOCATE(QETA) DEALLOCATE(UETA) DEALLOCATE(VETA) c l=index(init_out//' ',' ')-1 if (init_out(l:l) .ne. '/') then l=l+1 init_out(l:l)='/' endif write(fname,'(i7.7)') idat(5) open (1,file=init_out(1:l)//'preproc.bc.'//fname . ,status='unknown',form='unformatted') write(1) run,idat(1),idat(2),idat(3),idat(4),pdb,tb,qb,ub,vb write(6,*) 'in interp, pdb(1,1),pdb(1,2): ', pdb(1,1), + pdb(1,2) close(1) 671 format('L,U,V(LEV=20) ',I3,1x,f6.2,1x,f6.2,1x,f6.2) write(6,*) 'took square of neg number at ', negsum,' points' c c return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC subroutine levext(array,idim,jdim,kdim) real array(idim,jdim,kdim) do K=1,KDIM armax=-9999999. armin=9999999. imax=-9999. jmax=-9999. C do J=1,jdim C do I=1,idim do J=2,jdim-1 do I=2,idim-1 if (array(I,J,K).gt.armax) then armax=array(I,J,K) imax=i jmax=j endif if (array(I,J,K).lt.armin) armin=array(I,J,K) enddo enddo write(6,*) 'extremes for LEV: ', K, 'are ', armin,armax, + 'imax,jmax ', + imax,jmax enddo return end