program degrib2model_driver c implicit none c integer nx,ny,nz,kpds(200),datsav(9),n,k !Degrib grid dimensions c character(LEN=255):: gribfile,outfile c real esat,es c common /estab/esat(15000:45000),es(15000:45000) c_______________________________________________________________________________ c c *** Initialize tables. c call es_ini c c *** Read file names. c read(5,'(a)') gribfile read(5,'(a)') outfile c! gribfile='../data_in/init/gfsdata'//c_k//'.grb' c! outfile ='../data_in/init/gfsdata'//c_k//'.dat' c n=index(gribfile,' ')-1 print *,n write(6,*) 'first call of GET_GDS' c CALL GET_GDS(gribfile(1:n),datsav,kpds)!!!!!!!!!!!!!!!!!!! DRAGAN c nx=datsav(1) !!!!!!!!!!!!!!! c ny=datsav(2) !!!!!!!!!!!!!!!!!! nx=360 !!!!!!!!!!!!!!!!!! ny=181 !!!!!!!!!!!!!!!!!! if (nx .eq. 360 .and. ny .eq. 181) then nz=26 elseif (nx .eq. 144 .and. ny .eq. 73) then nz=17 else write(6,*) 'what grid is this??? ', nx,ny write(6,*) ' ' write(6,*) 'I do not know the vertical levels for this ' write(6,*) 'grid...please add some code to dgeta2model.f' write(6,*) ' ' STOP 999 endif print *,"nx,ny,nz,",nx,ny,nz,gribfile,outfile c call degrib2model(gribfile,outfile,nx,ny,nz) c end c c=============================================================================== c subroutine degrib2model(gribfile,outfile,nx,ny,nz) c implicit none c real cp,kappa,dlat,dlon,sum,avg parameter(cp=1004.,kappa=287./1004.) c integer nx,ny,nz,i,j,k,l,n,it,ip,jp,nsfcfld . ,iyear,imonth,iday,ifcsthr . ,count,kgds(200),ireto,n1,n2,n3,LenStr c real ht(nx,ny,nz) !Isobaric heights (m) . ,tp(nx,ny,nz) !Isobaric temps (K) . ,th(nx,ny,nz) !Isobaric theta (K) . ,uw(nx,ny,nz) !Isobaric u-wind (m/s) . ,vw(nx,ny,nz) !Isobaric v-wind (m/s) . ,rh(nx,ny,nz) !Isobaric rh,mr (%,kg/kg) . ,pr(nz),pri(nz) !Isobaric pressures (mb) . ,lat1,lat2,lon0,sw(2),ne(2) . ,xe,mrsat,esat,es . ,avnlevs(26),reanlevs(17) . ,slp(nx,ny,14) !slp(i,j,1)=EMSL;slp(i,j,2)=PMSL !@@@ . ,temp(nx,ny) !slp(i,j,3)=PSFC;slp(i,j,4)=ZSFC !slp(i,j,5 & 6) are 0-10 cm STC and SMC !slp(i,j,7 & 8) are 10-200 cm STC and SMC REAL, ALLOCATABLE:: htout(:,:,:), + uwout(:,:,:),vwout(:,:,:),rhout(:,:,:), + slpout(:,:,:) c character*255 gribfile,outdir,outfile,gdsfile character*2 gproj character*11 atime character*8 atime_r character*7 model character*5 timex c ! character*255 gribfile,outfile,gdsfile c ! character*2 gproj c ! character*11 atime c ! character*7 model data AVNLEVS/1000,975,950,925,900,850,800,750,700,650, + 600,550,500,450,400,350,300,250,200,150, + 100, 70, 50, 30, 20, 10/ data REANLEVS/1000,925,850,700,600,500, + 400,300,250,200,150,100, + 70, 50, 30, 20, 10/ c common /estab/esat(15000:45000),es(15000:45000) c_______________________________________________________________________________ c c *** Fill pressure levels. c Cmp nsfcfld will be 12 although 4 of these are set to missing for AVN... print *,'degrib2model...' nsfcfld=14 c nsfcfld=12 !@@@ Cmp write(6,*) 'using hardwired stuff! ', nx,ny,nz do k=1,nz if (nz .eq. 26) then pr(k)=AVNLEVS(k) elseif (nz .eq. 17) then pr(k)=reanlevs(k) endif C write(6,*) 'pr(K): ', pr(K) enddo c c *** Read in degrib data. c ccccccc n=index(gribfile,' ')-1 c call read_degrib(gribfile(1:n) c . ,nx*ny,nz,pr,ht,tp,rh,uw,vw,slp,atime)!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! write(6,*) 'return read_degrib' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ccc atime_r='05070400000' c open(12,file='./adate.txt') c read(12,*) atime_r c close(12) ccc n2=index(gribfile,' ')-1 ccc n1=index(gribfile,'.')+1 ccc n3=n2-n1+1 ccc LenStr=n3+8 ccc timex(1:n3)=gribfile(n1:n2) ccc write(6,*)n1,n2,LenStr,timex ccc atime(1:LenStr)=atime_r//timex(1:n3) ccc print *, 'ATIME =',ATIME ccc n=index(gribfile,' ')-1 cccccccccc write(*,*) 'OPEN FILE ',gribfile(1:n) ccccccccc open(12,file=gribfile(1:n),form='unformatted') ccccccc write(*,*) 'OPEN FILE ',gribfile(1:n) ccc OPEN(12,file='../../../2001070100.bin',form='unformatted', ccc & access='direct',status='old',recl=360*181*4) ccc do l=1,nz ccc read(12) ((ht(i,j,l),i=1,nx),j=ny,1,-1) ccc enddo ccc do l=1,nz ccc read(12) ((uw(i,j,l),i=1,nx),j=ny,1,-1) ccc enddo ccc do l=1,nz ccc read(12) ((vw(i,j,l),i=1,nx),j=ny,1,-1) ccc enddo ccc do l=1,nz ccc read(12) ((tp(i,j,l),i=1,nx),j=ny,1,-1) !!!!!!!!!! ccc enddo c do l=1,nz-1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ccc do l=1,nz ccc read(12) ((rh(i,j,l),i=1,nx),j=1,ny) ccc enddo c# do l=1,12 !@@@ c# read(12) ((slp(i,j,l),i=1,nx),j=1,ny) c print *,'SLP',l,slp(5,5,l) ccc read(12) ((slp(i,j,2),i=1,nx),j=1,ny) ccc print *,'SLP',2,slp(100,100,2) ccc read(12) ((slp(i,j,1),i=1,nx),j=1,ny) ccc print *,'SLP',1,slp(100,100,1) ccc read(12) ((slp(i,j,3),i=1,nx),j=1,ny) ccc print *,'SLP',3,slp(100,100,1) ccc read(12) ((slp(i,j,4),i=1,nx),j=1,ny) ccc print *,'SLP',4,slp(100,100,1) ccc do l=5,12 ccc read(12) ((slp(i,j,l),i=1,nx),j=ny,1,-1) ccc print *,'SLP',l,slp(100,100,l) ccc enddo ccc read(12) ((slp(i,j,13),i=1,nx),j=1,ny) ccc print *,'SLP',13,slp(100,100,1) ccc read(12) ((slp(i,j,14),i=1,nx),j=1,ny) ccc print *,'SLP',14,slp(100,100,1) c OPEN(11,file='../../../2001070100.bin',form='unformatted', c & access='direct',status='old',recl=360*181*4) cccccccc OPEN(11,file='../../../../1997010100.bin',form='unformatted', ccccccccc & access='direct',status='old',recl=360*181*4) cccccccccc OPEN(11,file='../../../1998020900.bin',form='unformatted', cccccccccc & access='direct',status='old',recl=360*181*4) ccccc OPEN(11,file='../../../../1996020900.bin',form='unformatted', ccccc & access='direct',status='old',recl=360*181*4) OPEN(11,file='../INIT_DATA/feb1996/1996020900.bin', & form='unformatted', & access='direct',status='old',recl=360*181*4) c OPEN(11,file='../../INIT_DATA/1997/1997010100.bin',form='unformatted', c & access='direct',status='old',recl=360*181*4) c do 20 i=1,26,1 c j=i c read(11,rec=j) uw(:,:,i) c 20 continue c print*,'uw=',uw(1,1,1) c do 30 i=1,26,1 c j=i+26 c read(11,rec=j) vw(:,:,i) c 30 continue c print*,'vw=',vw(1,1,1) c do 40 i=1,26,1 c j=i+52 c read(11,rec=j) ht(:,:,i) c 40 continue c print*,'ht=',ht(1,1,1) c do 50 i=1,26,1 c j=i+78 c read(11,rec=j) tp(:,:,i) c 50 continue c print*,'tp=',tp(1,1,1) c do 60 i=26,1,-1 c j=i+104 c read(11,rec=j) rh(:,:,i) c 60 continue c print*,'rh=',rh(1,1,1) c read(11,rec=131) slp(:,:,1) c read(11,rec=132) slp(:,:,2) c read(11,rec=133) slp(:,:,3) c read(11,rec=134) slp(:,:,4) c read(11,rec=135) slp(:,:,5) c read(11,rec=136) slp(:,:,6) c read(11,rec=137) slp(:,:,7) c read(11,rec=138) slp(:,:,8) c read(11,rec=139) slp(:,:,9) c read(11,rec=140) slp(:,:,10) c read(11,rec=141) slp(:,:,11) c read(11,rec=142) slp(:,:,12) c read(11,rec=143) slp(:,:,13) c read(11,rec=144) slp(:,:,14) c############################################# do 20 l=1,nz k=l read(11,rec=k) ((uw(i,j,l),i=1,nx),j=ny,1,-1) 20 continue c print*,'uw=',uw(1,1,1) do 30 l=1,nz k=l+26 read(11,rec=k) ((vw(i,j,l),i=1,nx),j=ny,1,-1) 30 continue c print*,'vw=',vw(1,1,1) do 40 l=1,nz k=l+52 read(11,rec=k) ((ht(i,j,l),i=1,nx),j=ny,1,-1) 40 continue c print*,'ht=',ht(1,1,1) do 50 l=1,nz k=l+78 read(11,rec=k) ((tp(i,j,l),i=1,nx),j=ny,1,-1) 50 continue c print*,'tp=',tp(1,1,1) do 60 l=1,nz k=l+104 read(11,rec=k) ((rh(i,j,l),i=1,nx),j=ny,1,-1) 60 continue c print*,'rh=',rh(1,1,1) read(11,rec=131) ((slp(i,j,1),i=1,nx),j=1,ny) read(11,rec=132) ((slp(i,j,2),i=1,nx),j=1,ny) read(11,rec=133) ((slp(i,j,3),i=1,nx),j=1,ny) read(11,rec=134) ((slp(i,j,4),i=1,nx),j=1,ny) read(11,rec=135) ((slp(i,j,5),i=1,nx),j=ny,1,-1) read(11,rec=136) ((slp(i,j,6),i=1,nx),j=ny,1,-1) read(11,rec=137) ((slp(i,j,7),i=1,nx),j=ny,1,-1) read(11,rec=138) ((slp(i,j,8),i=1,nx),j=ny,1,-1) read(11,rec=139) ((slp(i,j,9),i=1,nx),j=ny,1,-1) read(11,rec=140) ((slp(i,j,10),i=1,nx),j=ny,1,-1) read(11,rec=141) ((slp(i,j,11),i=1,nx),j=ny,1,-1) read(11,rec=142) ((slp(i,j,12),i=1,nx),j=ny,1,-1) ccc read(11,rec=143) ((slp(i,j,13),i=1,nx),j=ny,1,-1) ccc read(11,rec=144) ((slp(i,j,14),i=1,nx),j=ny,1,-1) OPEN(77,file='../INIT_DATA/feb1998/ICEC_WEASD.bin', & form='unformatted',access='direct',status='old',recl=360*181*4) read(77,rec=1) ((slp(i,j,13),i=1,nx),j=ny,1,-1) read(77,rec=2) ((slp(i,j,14),i=1,nx),j=ny,1,-1) do i=1,nx do j=1,ny if (slp(i,j,13).eq.9.9990003E+20) then slp(i,j,13)=0. endif if (slp(i,j,14).eq.9.9990003E+20) then slp(i,j,14)=0. endif enddo enddo do i=1,nx doj=1,ny ! print*,slp(i,j,13) enddo enddo c######################################### print *,'SLP',1,slp(234,70,1) print *,'SLP',2,slp(234,70,2) print *,'SLP',3,slp(234,70,3) print *,'SLP',4,slp(234,70,4) print *,'SLP',5,slp(234,70,5) print *,'SLP',6,slp(234,70,6) print *,'SLP',7,slp(234,70,7) print *,'SLP',8,slp(234,70,8) print *,'SLP',9,slp(234,70,9) print *,'SLP',10,slp(234,70,10) print *,'SLP',11,slp(234,70,11) print *,'SLP',12,slp(234,70,12) print *,'SLP',13,slp(234,70,13) print *,'SLP',14,slp(234,70,14) do l=5,12 c! print *,'SLP',14,slp(:,:,l) enddo c print*,slp(:,:,5) c print*,rh(22,:,23) close(11) close(77) ccc close(12) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c do k=1,nz c do i=1,nx c do j=1,ny c! print *,i,j,k,rh(i,j,k) c enddo c enddo c enddo print *,'rh(1,1,1),',rh(1,1,1) c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! do i=1,nx do j=ny,1,-1 do l=1,nz !dragan15 if (rh(i,j,l).eq.9.9990003E+20) print*,rh(i,j,l),i,j,l enddo do l=1,nz if (ht(i,j,l).eq.9.9990003E+20) ht(i,j,l)=-99999. if (uw(i,j,l).eq.9.9990003E+20) uw(i,j,l)=-99999. if (vw(i,j,l).eq.9.9990003E+20) vw(i,j,l)=-99999. if (tp(i,j,l).eq.9.9990003E+20) tp(i,j,l)=-99999. if (rh(i,j,l).eq.9.9990003E+20) rh(i,j,l)=-99999. enddo do l=1,14 !!!!!!!!!!!??? cccccc if (slp(i,j,l).eq.9.9990003E+20) slp(i,j,l)=0. enddo c############################################## DRAGAN do l=22,nz c if (rh(i,j,l).gt.10) rh(i,j,l)=10 !!!!!!!!!!!!!!!!????????????????????????????????? enddo do l=1,21 if (rh(i,j,l).gt.100) rh(i,j,l)=100 enddo do l=1,nz if (rh(i,j,l).lt.0) rh(i,j,l)=1 enddo c if(slp(i,j,1).ge.0.51) then c slp(i,j,1)=1 c else c slp(i,j,1)=0 c endif enddo enddo c########################################################################### c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c print*,rh(22,:,15) print*,slp(360,:,3) c c *** Check for any missing data. c do k=1,nz if (ht(1,1,k) .eq. -99999.) then print *,'Height data missing at level: ',pr(k) stop elseif (tp(1,1,k) .eq. -99999.) then print *,'Temperature data missing at level: ',pr(k) stop elseif (rh(1,1,k) .eq. -99999.) then print *,'RH data missing at level: ',k,pr(k) print *,'Calling RH patch.' call rh_fix(nx,ny,nz,rh,pr) elseif (uw(1,1,k) .eq. -99999.) then print *,'U-wind data missing at level: ',pr(k) stop elseif (vw(1,1,k) .eq. -99999.) then print *,'V-wind data missing at level: ',pr(k) stop endif enddo Cmp Cmp Handle missing surface data.... Cmp C C For reanalysis runs these data will not be used. Just bogus C to hold a place. C write(6,*) 'going to use nsfcfld= ', nsfcfld C If rean put bogus data into slp if (nz .eq. 17) then do K=5,8 do j=1,ny do i=1,nx if (mod(k,2) .eq. 0.and.k.ne.2) then !!!!!!!!!!!????????????????? slp(i,j,k)=0.14 else slp(i,j,k)=275. endif enddo enddo enddo endif !!! do k=5,nsfcfld do k=5,12 if (slp(1,1,k) .eq. -99999.) then write(6,*) 'SURFACE DATA MISSING... ', K if (k.eq.9.or.k.eq.11) then write(6,*) 'filling soil temp data...' do j=1,ny do i=1,nx slp(i,j,k)=slp(i,j,7) enddo enddo elseif(k.eq.10.or.k.eq.12.) then write(6,*) 'filling soil moisture data...' do j=1,ny do i=1,nx slp(i,j,k)=slp(i,j,8) enddo enddo endif endif Cmp if (mod(k,2) .eq. 0.and.k.ne.2) then !!!!!!!!!!!!!!!! do j=1,ny do i=1,nx if (slp(i,j,k) .eq. 0) then slp(i,j,k)=0.14 endif enddo enddo else do j=1,ny do i=1,nx if (slp(i,j,k) .eq. 0) then slp(i,j,k)=273.15 endif enddo enddo endif Cmp enddo c c *** Convert 3d temp to theta. c *** Compute Exner function. c *** Convert 3d rh to mr. c do k=1,nz pri(k)=1./pr(k) enddo c do k=1,nz do j=1,ny do i=1,nx th(i,j,k)=tp(i,j,k)*(1000.*pri(k))**kappa C ex(i,j,k)=cp*tp(i,j,k)/th(i,j,k) it=tp(i,j,k)*100 it=min(45000,max(15000,it)) xe=esat(it) mrsat=0.00622*xe/(pr(k)-xe) rh(i,j,k)=rh(i,j,k)*mrsat enddo enddo enddo Cmp Cmp for some reason AVN data appears to be "upside down" in the N-S Cmp sense. Flip the arrays of data. If reanalysis data, skip this part Cmp c! if (nz .eq. 17) goto 969 do k=1,nz do j=1,ny do i=1,nx temp(i,j)=th(i,j,k) enddo enddo do J=1,ny do i=1,nx th(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx temp(i,j)=uw(i,j,k) enddo enddo do J=1,ny do i=1,nx uw(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx temp(i,j)=vw(i,j,k) enddo enddo do J=1,ny do i=1,nx vw(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx temp(i,j)=ht(i,j,k) enddo enddo if (k.eq.1) write(6,*) 'flipping isobaric data - ' do J=1,ny do i=1,nx ht(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx temp(i,j)=rh(i,j,k) enddo enddo do J=1,ny do i=1,nx rh(i,j,k)=temp(i,ny-J+1) enddo enddo C====================================== do j=1,ny do i=1,nx C temp(i,j)=ex(i,j,k) enddo enddo do J=1,ny do i=1,nx C ex(i,j,k)=temp(i,ny-J+1) enddo enddo C----------------------------------------- enddo Cmp generalize eventually! ??????????????????????????????????????????????????? do k=1,nsfcfld C write(6,*) 'flipping surface fields... ' do j=1,ny do i=1,nx temp(i,j)=slp(i,j,k) enddo enddo do J=1,ny do i=1,nx slp(i,j,k)=temp(i,ny-J+1) enddo enddo enddo 969 continue CC CC CC Add bogus column if reanalysis data CC CC if (nx .eq. 144) then ALLOCATE(HTOUT(145,ny,nz),UWOUT(145,ny,nz),VWOUT(145,ny,nz)) ALLOCATE(RHOUT(145,ny,nz),SLPOUT(145,ny,14)) DO J=1,NY DO I=1,NX DO K=1,NSFCFLD SLPOUT(I,J,K)=SLP(I,J,K) ENDDO DO K=1,NZ HTOUT(I,J,K)=HT(I,J,K) UWOUT(I,J,K)=UW(I,J,K) VWOUT(I,J,K)=VW(I,J,K) RHOUT(I,J,K)=RH(I,J,K) ENDDO ENDDO DO K=1,NSFCFLD SLPOUT(145,J,K)=SLP(1,J,K) ENDDO DO K=1,NZ HTOUT(145,J,K)=HT(1,J,K) UWOUT(145,J,K)=UW(1,J,K) VWOUT(145,J,K)=VW(1,J,K) RHOUT(145,J,K)=RH(1,J,K) ENDDO ENDDO endif c c c *** Create output file name. c n=index(outfile,' ')-1 Cmp make ETA_avn to differentiate from files based on eta grib input if (nz .eq. 26) then model='ETA_avn' elseif (nz .eq. 17) then model='ETA_rnl' nx=145 endif n=index(outfile,' ')-1 cccccccccccccc outfile=outdir(1:n)//'/'//atime(1:LenStr)//'.'//model print *,model,' data ---> ',outfile(1:n) !@@@ open(1,file=outfile(1:n),status='unknown',form='unformatted') c c *** Write header stuff. c ip=1 jp=1 Cmp nsfcfld will be 12 although 4 of these are set to missing for AVN... c nsfcfld=12 !@@@ nsfcfld=14 gproj='LL' write(1) nz C new stuff for GDS file C C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C new stuff for GDS file C n=index(gribfile,' ')-1 write(6,*) 'calling get_fullgds with ', + gribfile(1:n) c call get_fullgds(gribfile(1:n),nx,ny,kgds) kgds(2)=360 kgds(3)=181 kgds(4)=-90000 kgds(6)=128 kgds(7)=90000 kgds(8)=-1000 kgds(9)=1000 kgds(10)=1000 kgds(20)=255 c!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! c n=index(gribfile,' ')-1 c write(6,*) 'calling get_fullgds with ', c + gribfile(1:n) !!!!!!commented DRAGAN c call get_fullgds(gribfile(1:n),nx,ny,kgds) write(6,*) 'back from get_fullgds' c call baclose(11,ireto) c if (nx .eq. 145) then c KGDS(8)=0 !@@@ ???? c KGDS(2)=145 c endif n=index(outfile,' ')-1 gdsfile='gdsinfo.ETA_avn' n=index(gdsfile,' ')-1 write(6,*) 'GDS written to ', gdsfile(1:n) open(unit=14,file=gdsfile(1:n),form='unformatted', + access='sequential') C C MODIFY GDS TO REFLECT FLIPPING OF DATA IN N-S SENSE C if (nz .eq. 26) then kgds(4)= -kgds(4) kgds(7)= -kgds(7) endif write(6,*) 'writing kgds ', (kgds(I),i=1,14) write(14) kgds close(14) ! write(6,*) 'past 14 close' C C end new stuff c do i=1,nx c do j=ny,1,-1 c if(slp(i,j,1).ne.0.and.slp(i,j,1).ne.1) then c if(slp(i,j,2).eq.0) then c print*,slp(i,j,1),i,j c endif C print*,slp(360,j,5) c enddo c enddo c print*,slp(:,:,3) c c *** Write isobaric upper air data. c write(6,*) 'starting writes to unit1' if (nx .ne. 145) then write(1) uw write(1) vw write(1) ht write(1) rh write(1) pr write(1) slp else WRITE(1) UWOUT WRITE(1) VWOUT WRITE(1) HTOUT WRITE(1) RHOUT WRITE(1) PR WRITE(1) SLPOUT endif if (allocated(uwout)) deallocate(uwout) if (allocated(vwout)) deallocate(vwout) if (allocated(htout)) deallocate(htout) if (allocated(rhout)) deallocate(rhout) if (allocated(slpout)) deallocate(slpout) Cmp write(6,*) 'writing these sfc values to the ETA_avn file...' do k=1,nsfcfld write(6,*) 'field, value ', k, slp(nx/2,ny/2,k) enddo do k=1,nsfcfld write(6,*) 'field, value ', k, slp(3*nx/4,3*ny/4,k) enddo c close(1) c return end c c=============================================================================== c subroutine rh_fix(nx,ny,nz,rh,pr) c implicit none c integer*4 nx,ny,nz,i,j,k,kk c real*4 rh(nx,ny,nz),pr(nz) c_______________________________________________________________________________ c c *** Fix bottom levels if necessary. c if (rh(1,1,1) .eq. -99999.) then do k=2,nz do j=1,ny do i=1,nx ! print *,i,j,k,rh(i,j,k) enddo enddo if (rh(1,1,k) .ne. -99999.) then do kk=k-1,1,-1 print *,'Copying',nint(pr(kk+1)),' mb to' . , nint(pr(kk)),' mb.' do j=1,ny do i=1,nx rh(i,j,kk)=rh(i,j,kk+1) enddo enddo enddo goto 10 endif enddo print *,'RH patch did not work.' stop endif c c *** Fix upper levels if necessary. c 10 continue if (rh(1,1,nz) .eq. -99999.) then do k=nz-1,1,-1 if (rh(1,1,k) .ne. -99999.) then do kk=k+1,nz print *,'Copying',nint(pr(kk-1)),' mb to' . , nint(pr(kk)),' mb.' do j=1,ny do i=1,nx rh(i,j,kk)=rh(i,j,kk-1) enddo enddo enddo goto 20 endif enddo endif c 20 continue do k=1,nz if (rh(1,1,k) .eq. -99999.) then print *,'RH patch did not work.' stop endif enddo c return end c c=============================================================================== c subroutine es_ini c common /estab/esat(15000:45000),es(15000:45000) c c *** Create tables of the saturation vapour pressure with up to c two decimal figures of accuraccy: c do it=15000,45000 t=it*0.01 p1 = 11.344-0.0303998*t p2 = 3.49149-1302.8844/t c1 = 23.832241-5.02808*alog10(t) esat(it) = 10.**(c1-1.3816E-7*10.**p1+ . 8.1328E-3*10.**p2-2949.076/t) es(it) = 610.78*exp(17.269*(t-273.16)/(t-35.86)) enddo c return end c c=============================================================================== CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC