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 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=alog(pt) 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 N=1,LM ETA(N+1)=ETA(N)+DETA(N) ENDDO eta(lmp1)=1.0 c c *** Compute ldt1. c dolm: do l=1,lm i=1 ii=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 do while (peta(l) .ge. spl(i) .and. i .lt. ldm) peta(l)=(prf0-pt)*eta(l)+pt i=i+1 enddo if (.not.lInterpAllLev) then do ii=1, ldm if ((abs(peta(l)-spl(ii)))*1E-2 .lt.2 )then MaskInt(l)=1+ii iii=ii exit endif MaskInt(l)=0 iii=i-1 enddo else MaskInt(l)=0 endif C if (i .le. 2) then ldt1(l)=1 elseif (i .eq. ldm) then ldt1(l)=ldm-2 else dpu=abs(spl(i-2)-peta(l)) dpd=abs(spl(i+1)-peta(l)) if (dpu .gt. dpd) then ldt1(l)=i-1 else ldt1(l)=i-2 endif endif enddo dolm c c *** Compute half eta values. c do l=1,lm etal(l)=0.5*(eta(l)+eta(l+1)) enddo 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 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 point's 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 point's 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) do ld=1,ldmm1 alpr(ld)=alog(spl(ld+1)/spl(ld)) enddo do ld=1,ldm Cmp ln of preset pressure levels alp(ld)=alog(spl(ld)) alpsq(ld)=alp(ld)**2 enddo do l=1,lm l1=ldt1(l) Cmakesworse dlt(l)=alpr(l1+1)*alog(spl(l1+2)/spl(l1))*(-alpr(l1)) dlt(l)=alpr(l1+1)*(alp(l1+2)-alp(l1))*(-alpr(l1)) enddo 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 l=min(lm,lhgt(i,j)) l1=ldt1(l) l2=l1+1 l3=l1+2 h1b=h(i,j,l1) h2 =h(i,j,l2) h3 =h(i,j,l3) h3m2=h3-h2 h3m1=h3-h1b b=(h3m2*(alpsq(l3)-alpsq(l1)) . -h3m1*(alpsq(l3)-alpsq(l2)))/dlt(l) c=(h3m1*(alp(l3)-alp(l2)) . -h3m2*(alp(l3)-alp(l1)))/dlt(l) if (abs(c) .lt. 0.1) then write(6,*) 'small abs(c).... ', i,j,c write(6,*) l,l1,l2,l3,h3m1,h3m2 write(6,*) alp(l1),alp(l3),alp(l2) write(6,*) dlt(l) write(6,*)h(i,j,l1), h(i,j,l2),h(i,j,l3) goto 532 endif Cmp if (b**2-4.*c*(-c*alpsq(l3)-b*alp(l3)+h3-hgt(i,j)).lt.0) then negsum=negsum+1 goto 532 endif alpgt=(-b-sqrt(b**2-4.*c . *(-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 pd(i,j)=(exp(alpgt)-pt)/ref(i,j) 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 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ write(6,*) 'special case.... ', i,j,pd(i,j) 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 Cmp if (i .ge. IM-1) write(6,*) 'i,j,pd,ref ', i,j,pd(i,j),ref(i,j) Cmp reinstate Cremove endif 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 dolm1: do l=lm,1,-1 doout2: do j=1,jm doin2: do i=1,im l1=ldt1(l) l2=l1+1 l3=l1+2 h1b=h(i,j,l1) h2 =h(i,j,l2) h3 =h(i,j,l3) h3m2=h3-h2 h3m1=h3-h1b b=(h3m2*(alpsq(l3)-alpsq(l1)) . -h3m1*(alpsq(l3)-alpsq(l2)))/dlt(l) c=(h3m1*(alp(l3)-alp(l2)) . -h3m2*(alp(l3)-alp(l1)))/dlt(l) c ************ Computation of heights at upper boundaries of eta layers. c q1 =q(i,j,l1) q2b=q(i,j,l2) q3 =q(i,j,l3) q3m2=q3-q2b q3m1=q3-q1 bq=(q3m2*(alpsq(l3)-alpsq(l1)) . -q3m1*(alpsq(l3)-alpsq(l2)))/dlt(l) cq=(q3m1*(alp(l3)-alp(l2)) . -q3m2*(alp(l3)-alp(l1)))/dlt(l) Chou inc Tp T1 =T(i,j,l1) T2b=T(i,j,l2) T3 =T(i,j,l3) T3m2=T3-T2b T3m1=T3-T1 bT=(T3m2*(alpsq(l3)-alpsq(l1)) . -T3m1*(alpsq(l3)-alpsq(l2)))/dlt(l) cT=(T3m1*(alp(l3)-alp(l2)) . -T3m2*(alp(l3)-alp(l1)))/dlt(l) Chou in Tp C alpeta=alog(pt+pd(i,j)*eta(l)) alpeta=alog(pt+pd(i,j)*eta(l)) alpq=alog(pt+pd(i,j)*(eta(l)+0.5*deta(l))) hetaij=h2+b*(alpeta-alp(l2))+c*(alpeta**2-alpsq(l2)) Ctest heta(i,j,l)=amax1(hetaij,0.01) if (lInterpAllLev) then heta(i,j,l) =hetaij CJLG Teste alteraçao Chou 202429 if (.not.lInterpTQLin) then qeta(i,j,l) = q2b+bq*(alpq-alp(l2)) . +cq*(alpq**2-alpsq(l2)) Teta(i,j,l) = T2b+bT*(alpq-alp(l2)) . +cT*(alpq**2-alpsq(l2)) !chou in Tp endif CJLG Teste alteraçao Chou 202429 else if (MaskInt(l).gt.0 ) then write(9500,*)"No Vert inter at level L=", . l,lInterpAllLev lin=MaskInt(l)-1 heta(i,j,l) = h(i,j,lin) CJLG Teste alteraçao Chou 202429 if (.not.lInterpTQLin) then qeta(i,j,l) = q(i,j,lin) Teta(i,j,l) = T(i,j,lin) endif CJLG Teste alteraçao Chou 202429 else heta(i,j,l) =hetaij CJLG Teste alteraçao Chou 202429 if (.not.lInterpTQLin) then qeta(i,j,l) = q2b+bq*(alpq-alp(l2)) . +cq*(alpq**2-alpsq(l2)) Teta(i,j,l) = T2b+bT*(alpq-alp(l2)) . +cT*(alpq**2-alpsq(l2)) !chou in Tp endif CJLG Teste alteraçao Chou 202429 endif endif enddo doin2 enddo doout2 enddo dolm1 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) C write(6,*) 'I,J,pd,ref ',i,j,pd(i,j),ref(i,j),pdmax 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 C do 519 j=1,jm C ivl=1 C ivh=im-mod(j,2) Cmytest ivh=im C do 519 i=ivl,ivh C if (j .eq. 1 .or. j .eq. jm) goto 520 Corig if ( (i .eq.1 .and. mod(j+1,2).eq.0) .or. i.eq.im) goto 521 C if ( (i .eq.1 .and. mod(j,2).eq.0) .or. i.eq.im) goto 521 C pdvp=0.25*(pd(i,j-1)+pd(i+ivw(j),j)+pd(i+ive(j),j)+pd(i,j+1)) C goto 522 C 520 pdvp = 0.50*(pd(i+ivw(j),j)+pd(i+ive(j),j)) C goto 522 C 521 pdvp = 0.50*(pd(i,j-1)+pd(i,j+1)) 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' DO J=1,JM DO I=1,IM NN=(J-1)*IM+I IDXX(NN)=I JDXX(NN)=J ENDDO ENDDO NUMK=IM*JM C DO NN=1,NK C IDX(NN)=IDXX(NN) C JDX(NN)=JDXX(NN) C ENDDO MNKNT=1 do400: DO L=1,LM C KNTX=MNKNT C KNTIN=KNTX C KP1=MIN(LM,KNTIN+1) C NUMK=0 C ALPETAI(I,J)=ALOG(PT+PDVP(I,J)*ETAL(L)) C cdir$ ivdep do524out: DO N=1,NUMK I=IDXX(N) J=JDXX(N) ALPETAI(I,J)=ALOG(PT+PDVP(I,J)*ETAL(L)) ALPETK=ALPETAI(I,J) do524in: do ld=2,ldm if (lInterpAllLev) then if (ALPETK .GT. ALP(LD) .AND. LD .LT. LDM) . CYCLE do524in ULD=U(I,J,LD) VLD=V(I,J,LD) CF=(ALP(LD)-ALPETK)/ALPR(LD-1) UKB=ULD+(U(I,J,LD-1)-ULD)*CF VKB=VLD+(V(I,J,LD-1)-VLD)*CF UETA(I,J,L)=UKB VETA(I,J,L)=VKB CJLG Teste alteraçao Chou 202429 if (lInterpTQLin) then !Chou 20240929 including Teta and Qeta to follow the same interpolation as the !winds (mind v and h points at different position's) TTLD=T(I,J,LD) QQLD=Q(I,J,LD) TKB=TTLD+(T(I,J,LD-1)-TTLD)*CF QKB=QQLD+(Q(I,J,LD-1)-QQLD)*CF TETA(I,J,L)=TKB QETA(I,J,L)=QKB !Chou 20240929 end of inclusion Teta and Qeta (remember to declare: TTLD, QQLD, endif CJLG Teste alteraçao Chou 202429 EXIT do524in else if (MaskInt(l).gt.0 ) then lin=MaskInt(l)-1 UETA(I,J,L)=U(I,J,lin) VETA(I,J,L)=V(I,J,lin) CJLG Teste alteraçao Chou 202429 if (lInterpTQLin) then TETA(I,J,L)=T(I,J,LD) QETA(I,J,L)=Q(I,J,LD) endif CJLG Teste alteraçao Chou 202429 EXIT do524in else if (ALPETK .GT. ALP(LD) .AND. LD .LT. LDM) . CYCLE do524in ULD=U(I,J,LD) VLD=V(I,J,LD) CF=(ALP(LD)-ALPETK)/ALPR(LD-1) UKB=ULD+(U(I,J,LD-1)-ULD)*CF VKB=VLD+(V(I,J,LD-1)-VLD)*CF UETA(I,J,L)=UKB VETA(I,J,L)=VKB CJLG Teste alteraçao Chou 202429 if (lInterpTQLin) then TTLD=T(I,J,LD) QQLD=Q(I,J,LD) TKB=TTLD+(T(I,J,LD-1)-TTLD)*CF QKB=QQLD+(Q(I,J,LD-1)-QQLD)*CF TETA(I,J,L)=TKB QETA(I,J,L)=QKB !Chou 20240929 end of inclusion Teta and Qeta (remember to declare: TTLD, QQLD, endif CJLG Teste alteraçao Chou 202429 EXIT do524in endif endif enddo do524in ENDDO do524out enddo do400 DEALLOCATE(PDVP) write(6,*) 'past PDVP deallocate' C do j=1,jm C ivh=im-mod(j,2) C do i=1,ivh C do l=1,lm C if (l .eq. 1) alplb=alpt C alpub=alplb C alplb=alog(pt+pdvp(i,j)*eta(l+1)) C alpeta=sqrt(0.5*(alpub**2+alplb**2)) C do 524 ld=2,ldm C if (alpeta .gt. alp(ld) .and. ld .lt. ldm) goto 524 C uld=u(i,j,ld) C vld=v(i,j,ld) C ldm1=ld-1 C cf=(alp(ld)-alpeta)/alpr(ldm1) C uetaij(l)=uld+(u(i,j,ldm1)-uld)*cf C vetaij(l)=vld+(v(i,j,ldm1)-vld)*cf C if (abs(uetaij(l)) .gt. 100.) then C write(6,*) uld, u(i,j,ldm1), cf C endif C goto 523 C 524 continue C 523 continue C do l=1,lm C ueta(i,j,l)=uetaij(l) C veta(i,j,l)=vetaij(l) C enddo C enddo C enddo c c 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 do ivi=1,lm l=lmp1-ivi if (lGeo2Temp) then phlb=phub phub=g*heta(i,j,l) teta(i,j,l)=-(phlb-phub)*(pt+etal(l)*pd(i,j)) + /(r*deta(l)*pd(i,j)) endif if (ivi .eq. lm) then tsavetmp=teta(i,j,lm) endif Cmp teta(i,j,l)=amin1(teta(i,j,l),325.) teta(i,j,l)=amax1(teta(i,j,l),150.) Cmp if (eta(l+1) .le. ref(i,j)) then !Chou 20240929 if (teta(i,j,l). gt. a3 ) then ! Teste above 0oC rseta=pq0/(pd(i,j)*etal(l)+pt) . *exp(a2*(teta(i,j,l)-a3)/(teta(i,j,l)-a4)) corr= 1./(1.- vtmpc1 * rseta) qseta= rseta * corr !chou 20240929 from mixing ratio to specific humidity else rseta=pq0/(pd(i,j)*etal(l)+pt) . *exp(ai2*(teta(i,j,l)-a3)/(teta(i,j,l)-ai4)) corr= 1./(1.- vtmpc2 * rseta) qseta= rseta * corr !chou 20240929 from mixing ratio to specific humidity endif !Chou 202409 rheta=qeta(i,j,l)/qseta endif rheta=min(rheta,tresh) qeta(i,j,l)=rheta*qseta 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 c C************************ c ************ Now redefine pd to have it equal to ps-pt. C************************** c pd(i,j)=ref(i,j)*pd(i,j) C Cmp write(6,301) pd(i,j) 301 format('redefined pd',f9.2) C**************************** C************************** enddo enddo DEALLOCATE(HETA) DEALLOCATE(HGT) write(6,*) 'heta deallocated' c c *** Set u, v and t equal to zero at points below the ground. c do j=2,jm-1 ihl=1+mod(j,2) ihh=im-1 ! Cmp i indices changed from zj code do i=ihl,ihh do l=1,lm Cmp Cmp---change this part???????????????????????????? if (eta(l+1) .gt. ref(i,j)) then csd teta(i,j,l)=0. csd qeta(i,j,l)=0. C teta(i,j,l)=0. C qeta(i,j,l)=0. ueta(i+ihw(j),j,l)=0. veta(i+ihw(j),j,l)=0. ueta(i+ihe(j),j,l)=0. veta(i+ihe(J),j,l)=0. ueta(i,j-1,l)=0. veta(i,j-1,l)=0. ueta(i,j+1,l)=0. veta(i,j+1,l)=0. endif C if(abs(ueta(i,j,l)).gt.100.) write(6,*) 'ueta ',i,j,l,ueta(i,j,l) C if(abs(veta(i,j,l)).gt.100.) write(6,*) 'veta ',i,j,l,veta(i,j,l) enddo enddo enddo c c *** The separation of the boundary values. c C C C write(6,*) 'sample wind values, lev=20 ' C write(6,*) 'U,V,spd(1,1) ', ueta(1,1,20),veta(1,1,20), C + (ueta(1,1,20)**2.+veta(1,1,20)**2.)**0.5 C write(6,*) 'U,V,spd(35,35) ', ueta(35,35,20),veta(35,35,20), C + (ueta(35,35,20)**2.+veta(35,35,20)**2.)**0.5 C write(6,*) 'U,V,spd(im,4) ', ueta(im,4,20),veta(im,4,20), C + (ueta(im,4,20)**2.+veta(im,4,20)**2.)**0.5 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' C write(6,*) 'temps' C call levext(teta,im,jm,lm) C write(6,*) 'uwnd' C call levext(ueta,im,jm,lm) C write(6,*) 'vwnd' C call levext(veta,im,jm,lm) C write(6,*) 'Z' C call levext(heta,im,jm,lm) C write(6,*) 'pd' C call levext(pd,im,jm,1) C write(6,*) 'pd' do J=jm,jm/2,-1 C write(6,666) (pd(i,j)/100.,i=im-10,im) enddo C write(6,*) 'ueta at lm=19' do J=jm,jm/2,-1 C write(6,666) (ueta(i,j,19),i=im-10,im) enddo C write(6,*) 'veta at lm=19' do J=jm,jm/2,-1 C write(6,666) (veta(i,j,19),i=im-10,im) enddo C write(6,*) 'HMASK (NE corner) ' do J=jm,jm/2,-1 C write(6,667) (hmask(i,j),i=im-10,im) enddo C write(6,*) 'VMASK (NE corner) ' do J=jm,jm/2,-1 C write(6,667) (vmask(i,j),i=im-10,im) enddo 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 Ctst C open(unit=4,file=init_out(1:l)//'eta_temps.fil', C + status='unknown',form='unformatted') C write(4) teta Ctst write(6,*) 'ueta, veta, level 20' 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 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