!############################# Change Log ################################## ! 5.0.0 ! !########################################################################### ! Copyright (C) 1990, 1995, 1999, 2000, 2003 - All Rights Reserved ! Regional Atmospheric Modeling System - RAMS !########################################################################### subroutine chem_pressure_stage(n1,n2,nhem,glat,glon,glat2,glon2) use isan_coms !--(DMK-CCATT-INI)---------------------------------------------------------------- !srf-chem use chem_isan_coms use chem1_list, only: chemical_mechanism ! intent(in) use mem_chem1, only: CHEM_ASSIM, CHEMISTRY ! intent(in) !srf-chem-end !--(DMK-CCATT-FIM)---------------------------------------------------------------- use rconstants implicit none integer :: n1,n2,nhem real :: glat(n1,n2),glon(n1,n2),glat2(n1,n2),glon2(n1,n2) real :: fnprx,fnpry,grx,gry,gglat,gglon,thmax,thmin & ,xswlon_east,cntlon_east,rr,each_levp_real(maxpr) integer :: i,j,k,lv,ifm,loop,n,iunit real, external :: rs !--(DMK-CCATT-INI)---------------------------------------------------------------- !srf-chem character(len=32) :: chemical_mechanism_test !srf-chem-end !--(DMK-CCATT-FIM)---------------------------------------------------------------- ! Read the header of input pressure file. print 91,innpr(1:len_trim(innpr)) 91 format(//,' Reading pressure gridded data',/,a,//) !write(*,fmt='(2(A3,1X),4(A7,1X))') 'n1','n2','glat','glon','glat2','glon2' !write(*,fmt='(2(I3.3,1X),4(F7.2,1X))') n1,n2,glat(n1,n2),glon(n1,n2),glat2(n1,n2),glon2(n1,n2) open(11,file=innpr(1:len_trim(innpr))) read(11,*) marker,isversion if(marker.ne.999999) isversion=1 ! print *,'File to be read:',innpr(1:len_trim(innpr)),isversion select case (isversion) case (1) rewind 11 read(11,*) iyy,imm,idd,ihh,nprz,nprx,npry,xswlon,xswlat,gdatdx,gdatdy read(11,*) (each_levp_real(n),n=1,nprz) !dmae - geraDp generating integer, dprep generating real do n=1,nprz levpr(n) = int(each_levp_real(n)) end do write(*,*) "Info 1: ",iyy,imm,idd,ihh,nprz,nprx,npry,xswlon,xswlat,gdatdx,gdatdy write(*,*) "Levels: ",(levpr(n),n=1,nprz) inproj=1 if(iyy.lt.100) iyy=iyy+1900 ihh=ihh*100 case (2) print*,'doing RALPH 2 format' read(11,*) iyy,imm,idd,ihh,itinc,nprz,nprx,npry read(11,*) inproj,gdatdx,gdatdy,xswlat,xswlon & ,xnelat,xnelon,cntlat,cntlon,secondlat read(11,*) ivertcoord,(levpr(lv),lv=1,nprz) end select !--(DMK-CCATT-INI)---------------------------------------------------------------- !srf-chem if(CHEM_ASSIM==1 .and. CHEMISTRY >= 0) then ! read chemical mechanism that dp file is prepared for read(11,*) chemical_mechanism_test ! check consistency if(trim( chemical_mechanism_test ) /= trim(chemical_mechanism)) then print*,'chemical mechanism not the same as expected' print*,'expected=',trim(chemical_mechanism(1:len_trim(chemical_mechanism))) print*,'read =',trim(chemical_mechanism_test(1:len_trim(chemical_mechanism_test))) stop 'wrong chem mechanism at chem_astp' else print*,'ISAN-CHEM data assimilation for chemical mechanism: ',& trim(chemical_mechanism(1:len_trim(chemical_mechanism))) endif endif !srf-chem-end !--(DMK-CCATT-FIM)---------------------------------------------------------------- ! Check for consistency between file parameters and namelist parameters if(iyy.ne.iyear.or.imm.ne.imonth & .or.idd.ne.idate.or.ihh.ne.ihour) then print*,'Pressure file dates not the same as namelist!' print*,'Year :',iyy,iyear print*,'Month:',imm,imonth print*,'Day :',idd,idate print*,'Hour :',ihh,ihour stop 'pr_dates' endif ! Check pressure data domain size and location if (inproj/=1.and.nhem>0) then print*,'You must input a lat-lon pressure grid ' & ,'to run a global simulation !!!!' stop 'glob-no-press' endif if (inproj.eq.1) then ! If necessary, convert longitude specification to values in the range ! [-180.,180.]. xswlon=mod(xswlon+900.,360.)-180. xnelon=mod(xnelon-900.,-360.)+180. if(xswlon.lt.-180..or.xswlon.ge.180. .or. & xnelon.lt.-180..or.xnelon.gt.180.01.or. & xswlat.lt.-90. .or.xnelat.gt.90.01) then print*,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' print*,'!!! BAD DOMAIN SPECIFICATION !!' print*,'!!! xswlat,xswlon - ',xswlat,xswlon print*,'!!! xnelat,xnelon-xswlon - ',xnelat,xnelon-xswlon print*,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' stop 'bad_domain' endif fnprx=float(nprx) fnpry=float(npry) ! Set "global domain" flags to determine whether 2 extra rows will be ! added to borders of the input gridded data. iglobew = 1 if doing ! full 360 degrees of longitude. iglobs = 1 if including the south ! pole and if iglobew = 1, and iglobn = 1 if including ! the north pole and if iglobew = 1. iglobew=0 iglobn=0 iglobs=0 if(fnprx*gdatdx.gt.359.9) then if(abs(xswlon+180.).gt.1.e-3) then print*,'When extracting all longitudes of pressure data,' print*,'xswlon must be set to -180.' !stop 'xswlon' endif iglobew=1 endif if(iglobew.eq.1.and.xswlat.lt.-89.9) iglobs=1 if(iglobew.eq.1.and.xnelat.gt. 89.9) iglobn=1 idatelin=0 if(xnelon.lt.xswlon+.01) idatelin=1 ! Search coarse grid points for any that are outside the bounds of ! the pressure data. If any are found, stop. loop=1 if(nhem>0) loop=2 do ifm = 1,loop do j = 1,n2 do i = 1,n1 if(ifm==1) then gglat=glat(i,j) gglon=glon(i,j) elseif(ifm==2) then gglat=glat2(i,j) gglon=glon2(i,j) endif gry = (gglat - xswlat) / gdatdy + 1. if(gry.lt.2-iglobs.or.gry.gt.fnpry+iglobn-1) then print*,'Model grid point latitude must be' print*,'at least 1 gridpoint inside' print*,'pressure data area' print*,'ifm,i,j,glat,xswlat,xnelat' & ,ifm,i,j,gglat,xswlat,xnelat stop 'isan 1: outside p-data 1' endif if(idatelin.eq.1.and.gglon.lt.xswlon) gglon=gglon+360. grx=(gglon-xswlon)/gdatdx+1. if(grx.lt.2-iglobew.or.grx.gt.fnprx+iglobew-1) then print*,'Model grid point longitude must be' print*,'at least 1 gridpoint inside' print*,'pressure data area' print*,'ifm,i,j,glon,xswlon,xnelon' & ,ifm,i,j,gglon,xswlon,xnelon,iglobew,iglobs,iglobn stop 'isan 1: outside p-data 2' endif enddo enddo enddo endif ! Deallocate memory for the pressure data if(allocated(p_u)) deallocate(p_u) if(allocated(p_v)) deallocate(p_v) if(allocated(p_t)) deallocate(p_t) if(allocated(p_z)) deallocate(p_z) if(allocated(p_r)) deallocate(p_r) if(allocated(p_ur)) deallocate(p_ur) if(allocated(p_vr)) deallocate(p_vr) if(allocated(p_lat)) deallocate(p_lat) if(allocated(p_lon)) deallocate(p_lon) if(allocated(p_slp)) deallocate(p_slp) if(allocated(p_sfp)) deallocate(p_sfp) if(allocated(p_sft)) deallocate(p_sft) if(allocated(p_snow)) deallocate(p_snow) if(allocated(p_sst)) deallocate(p_sst) !--(DMK-CCATT-INI)---------------------------------------------------------------- !srf-chem if(allocated(p_sc)) deallocate(p_sc) !srf-chem-end !--(DMK-CCATT-FIM)---------------------------------------------------------------- ! Allocate memory for the pressure data allocate(p_u(nprx,npry,nprz)) allocate(p_v(nprx,npry,nprz)) allocate(p_t(nprx,npry,nprz)) allocate(p_z(nprx,npry,nprz)) allocate(p_r(nprx,npry,nprz)) !--(DMK-CCATT-INI)---------------------------------------------------------------- !srf-chem if(CHEM_ASSIM == 1 .and. nspecies>0 ) allocate(p_sc(nprx,npry,nprz,nspecies)) !srf-chem-end !--(DMK-CCATT-FIM)---------------------------------------------------------------- ! p_ur,p_vr arrays for rotated winds used later allocate(p_ur(nprx,npry,nprz)) allocate(p_vr(nprx,npry,nprz)) allocate(p_lat(nprx,npry)) allocate(p_lon(nprx,npry)) allocate(p_slp(nprx,npry)) allocate(p_sfp(nprx,npry)) allocate(p_sft(nprx,npry)) allocate(p_snow(nprx,npry)) allocate(p_sst(nprx,npry)) ! Fill with missing in case they are not present p_slp (1:nprx,1:npry)=1.e30 p_sfp (1:nprx,1:npry)=1.e30 p_sft (1:nprx,1:npry)=1.e30 p_snow(1:nprx,1:npry)=1.e30 p_sst (1:nprx,1:npry)=1.e30 iunit=11 print 92 92 format(////,1x,70('*')/ & ,' Access coarse resolution pressure data',/,1X,70('*')///) do k=1,nprz pnpr(k)=levpr(k)*100. enddo ! Call routine to fill pressure arrays from the chosen dataset. !--(DMK-CCATT-INI)---------------------------------------------------------------- call chem_get_press (iunit) !--(DMK-CCATT-OLD)---------------------------------------------------------------- ! call get_press (iunit) !--(DMK-CCATT-FIM)---------------------------------------------------------------- !!!!!!!! Be careful !!!!!!!!! ! Check input humidity variable. Assume that if the max of the field is greater ! than 1.1 (allow for some machine roundoff), ! it is specific humidity (in g/kg) which needs to be converted to rh !!$if(maxval(p_r(1:nprx,1:npry,1:nprz)) > 1.1)then !!$ print*,'------------------------------------------------' !!$ print*,' Converting specific humidity to rh',nprx,npry,nprz !!$ print*,'------------------------------------------------' !!$ do k=1,nprz !!$ do j=1,npry !!$ do i=1,nprx !!$ rr=p_r(i,j,k)*.001/(1.-p_r(i,j,k)*.001) !!$ p_r(i,j,k)=rr/rs(pnpr(k),p_t(i,j,k)) !!$ if(p_r(i,j,k) <= 0.1) p_r(i,j,k)=.1 !!$! if(i.eq.1.and.j.eq.1)print *,k,tn(i,j,k),pnpr(k),rr,rn(i,j,k),zn(i,j,k) !!$ enddo !!$ enddo !!$ enddo !!$endif ! Find max-min theta at bottom and top levels thmax=1. thmin=1000. do j=1,npry do i=1,nprx if(p_t(i,j,nprz).lt.1.e20) & thmax=max(thmax,p_t(i,j,nprz)*(p00/pnpr(nprz))**rocp) if(p_t(i,j,1).lt.1.e20) & thmin=min(thmin,p_t(i,j,1)*(p00/pnpr(1))**rocp) enddo enddo print 300,levpr(1),thmin,levpr(nprz),thmax 300 format(//,' Minimum THETA at ',I4,' mb - ',F8.2/ & ,' Maximum THETA at ',I4,' mb - ',F8.2) ! Compute lat-lon at input pressure data points select case (inproj) case (1) do j=1,npry do i=1,nprx p_lat(i,j)=xswlat+(j-1)*gdatdy p_lon(i,j)=xswlon+(i-1)*gdatdx enddo enddo case (2) print*,'lamb-con:',xswlat,xswlon,cntlat,cntlon,gdatdx,gdatdy do j=1,npry do i=1,nprx call lc2ll(cntlat,cntlon,p_lat(i,j),p_lon(i,j) & ,float(i),float(j),xswlat,xswlon,gdatdx) ! print*,i,j,p_lat(i,j),p_lon(i,j) enddo enddo case (3) print*,'polar:',cntlat,cntlon,gdatdx,gdatdy,xswlat,xswlon xswlon_east=xswlon if(xswlon<0.) xswlon_east=xswlon+360. cntlon_east=cntlon if(cntlon<0.) cntlon_east=cntlon+360. do j=1,npry do i=1,nprx call w3fb07(float(i),float(j),xswlat,xswlon_east,gdatdx,cntlon_east & ,p_lat(i,j),p_lon(i,j)) ! print*,i,j,p_lat(i,j),p_lon(i,j) enddo enddo end select close(iunit) return !--(DMK-CCATT-INI)---------------------------------------------------------------- end subroutine chem_pressure_stage !--(DMK-CCATT-OLD)---------------------------------------------------------------- !end subroutine pressure_stage !--(DMK-CCATT-FIM)---------------------------------------------------------------- subroutine chem_pressure_stage_grib2(n1,n2,nhem,glat,glon,glat2,glon2) !--(DMK-CCATT-OLD)---------------------------------------------------------------- !subroutine pressure_stage(n1,n2,nhem,glat,glon,glat2,glon2) !--(DMK-CCATT-FIM)---------------------------------------------------------------- use isan_coms use ModDateUtils, only: & date_add_to !--(DMK-CCATT-INI)---------------------------------------------------------------- !srf-chem use chem_isan_coms use chem1_list, only: chemical_mechanism ! intent(in) use mem_chem1, only: CHEM_ASSIM, CHEMISTRY ! intent(in) !srf-chem-end !--(DMK-CCATT-FIM)---------------------------------------------------------------- use rconstants !use dump, only: & ! dumpMessage #ifdef GRIB2 use wgrib2api #endif implicit none include "constants.f90" integer :: n1,n2,nhem real :: glat(n1,n2),glon(n1,n2),glat2(n1,n2),glon2(n1,n2) real :: fnprx,fnpry,grx,gry,gglat,gglon,thmax,thmin & ,xswlon_east,cntlon_east,rr,each_levp_real(maxpr) integer :: i,j,k,lv,ifm,loop,n,iunit,stat real, external :: rs character (len=200) :: metadata character(len=10) :: varName(5) !# composed name of levels (:press:) character (len=300) :: grid_info character (len=99) :: invline real, allocatable :: var(:,:),lat(:,:),lon(:,:) integer :: nx,ny,thisHour integer :: iyyl,imml,iddl,ihhl character(len=3) :: cnz !--(DMK-CCATT-INI)---------------------------------------------------------------- !srf-chem character(len=32) :: chemical_mechanism_test !srf-chem-end !--(DMK-CCATT-FIM)---------------------------------------------------------------- call str2int(innpr(len_trim(innpr)-2:len_trim(innpr)),thisHour,stat) #ifdef GRIB2 ! iErrNumber = nf90_open(path = trim(innpr), mode = nf90_nowrite, ncid = ncid) ! if (status /= nf90_noerr) iErrNumber=dumpMessage(c_tty,c_yes,header,modelVersion & ! ,c_fatal,trim(innpr)//' not found. Please, check it!') ! Read the header of input pressure file. varName(1)=':'//trim(temperature_varname)//':' print 91,innpr(1:len_trim(innpr)) 91 format(//,' Reading pressure gridded netCDF data',/,a,//) !print 92,innpr(1:len_trim(innpr))//'.inv' !92 format(//,' using the inventory GRIB2 data',/,a,//) !write(*,fmt='(2(A3,1X),4(A7,1X))') 'n1','n2','glat','glon','glat2','glon2' !write(*,fmt='(2(I3.3,1X),4(F7.2,1X))') n1,n2,glat(n1,n2),glon(n1,n2),glat2(n1,n2),glon2(n1,n2) !Making inventory using a wgrib2 function grb2_mk_inv iErrNumber = grb2_mk_inv(innpr(1:len_trim(innpr)),innpr(1:len_trim(innpr))//'.inv') ! get number of first var levels using grb2_inq from wgrib2 lib nprz_grib2 = grb2_inq(innpr(1:len_trim(innpr)),innpr(1:len_trim(innpr))//'.inv' & ,trim(varname(1)),' mb:') allocate(slevs(nprz_grib2)) !If the total of levels ir greater than defined in RAMSIN use the defined if(nprz_grib2>z_max_level) then nprz=z_max_level else nprz=nprz_grib2 endif do i = 1,nprz_grib2 ! get pressure leves and levels labels using grb2_inq from wgrib2 lib iErrNumber=grb2_inq(trim(innpr),trim(innpr)//'.inv',trim(varname(1)),' mb:' & ,sequential=i-1,desc=metadata) if (iErrNumber.ne.1) stop 2 j = index(metadata,trim(varname(1))) + len(trim(varname(1))) k = index(metadata," mb:") + len(" mb:")-1 read(metadata(j:),*) levpr_grib2(i) slevs(i) = metadata(j-1:k) !print *,'metadata(j:): ',metadata(j:),slevs(i),levpr_grib2(i) enddo !Getting information in GRIB2 file. Using the first variable (U) iErrNumber = grb2_inq(trim(innpr),trim(innpr)//'.inv',trim(varName(1)) & , data2=var & , lat=lat, lon=lon, grid_desc=grid_info, desc=invline) !Closing grib2 file to save memory iErrNumber = grb2_free_file(trim(innpr)) nxGrib=size(var,1) nyGrib=size(var,2) allocate(mask(nxGrib,nyGrib)) ! Increments lat and lon gdatdy=(lat(1,nyGrib)-lat(1,1))/(nyGrib-1) gdatdx=(lon(nxGrib,1)-lon(1,1))/(nxGrib-1) ! !Make number of points lon and lat accordingly setup in RAMSIN nprx=((final_longitude-initial_longitude)/gdatdx)+1 npry=((final_latitude-initial_latitude)/gdatdy)+1 write(cnz,fmt='(I3.3)') nprz_grib2 write(*,fmt='("Inventory(1): ",3(I5,1X))') nxGrib,nyGrib,nprz_grib2 write(*,fmt='("Inventory(2): ",4(E18.4,1X))') lon(1,1),gdatdx,lat(1,1),gdatdy write(*,fmt='("Inventory(3): "'//cnz//'(F6.1,1X))') (levpr_grib2(i),i=1,nprz_grib2) ! !convert date string got from grib2 file to integer call str2int(invline(3:6),iyyl,stat) call str2int(invline(7:8),imml,stat) call str2int(invline(9:10),iddl,stat) call str2int(invline(11:12),ihhl,stat) !Setting initial lon and lat xswlon=initial_longitude-360.0 xswlat=initial_latitude ! inproj=1 !if(iyy.lt.100) iyy=iyy+1900 !ihh=(ihh+thisHour)*100 !print *,iyyl,imml,iddl,ihhl,real(thisHour*3600) call date_add_to(iyyl,imml,iddl,ihhl,real(thisHour*3600),'s' & ,iyy,imm,idd,ihh) ihh=ihh/100 ! ! write(*,*) "Info 1: ",iyy,imm,idd,ihh,nprz,nprx,npry,xswlon,xswlat,gdatdx,gdatdy ! write(*,*) "Levels: ",(levpr(n),n=1,nprz) ! Check for consistency between file parameters and namelist parameters if(iyy.ne.iyear.or.imm.ne.imonth & .or.idd.ne.idate.or.ihh.ne.ihour) then print*,'Pressure file dates not the same as namelist!' print*,'Year :',iyy,iyear print*,'Month:',imm,imonth print*,'Day :',idd,idate print*,'Hour :',ihh,ihour stop 'pr_dates' endif ! Check pressure data domain size and location if (inproj/=1.and.nhem>0) then print*,'You must input a lat-lon pressure grid ' & ,'to run a global simulation !!!!' stop 'glob-no-press' endif if (inproj.eq.1) then ! If necessary, convert longitude specification to values in the range ! [-180.,180.]. xswlon=mod(xswlon+900.,360.)-180. xnelon=mod(xnelon-900.,-360.)+180. if(xswlon.lt.-180..or.xswlon.ge.180. .or. & xnelon.lt.-180..or.xnelon.gt.180.01.or. & xswlat.lt.-90. .or.xnelat.gt.90.01) then print*,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' print*,'!!! BAD DOMAIN SPECIFICATION !!' print*,'!!! xswlat,xswlon - ',xswlat,xswlon print*,'!!! xnelat,xnelon-xswlon - ',xnelat,xnelon-xswlon print*,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' stop 'bad_domain' endif fnprx=float(nprx) fnpry=float(npry) ! Set "global domain" flags to determine whether 2 extra rows will be ! added to borders of the input gridded data. iglobew = 1 if doing ! full 360 degrees of longitude. iglobs = 1 if including the south ! pole and if iglobew = 1, and iglobn = 1 if including ! the north pole and if iglobew = 1. iglobew=0 iglobn=0 iglobs=0 if(fnprx*gdatdx.gt.359.9) then if(abs(xswlon+180.).gt.1.e-3) then print*,'When extracting all longitudes of pressure data,' print*,'xswlon must be set to -180.' !stop 'xswlon' endif iglobew=1 endif if(iglobew.eq.1.and.xswlat.lt.-89.9) iglobs=1 if(iglobew.eq.1.and.xnelat.gt. 89.9) iglobn=1 idatelin=0 if(xnelon.lt.xswlon+.01) idatelin=1 ! Search coarse grid points for any that are outside the bounds of ! the pressure data. If any are found, stop. loop=1 if(nhem>0) loop=2 do ifm = 1,loop do j = 1,n2 do i = 1,n1 if(ifm==1) then gglat=glat(i,j) gglon=glon(i,j) elseif(ifm==2) then gglat=glat2(i,j) gglon=glon2(i,j) endif gry = (gglat - xswlat) / gdatdy + 1. if(gry.lt.2-iglobs.or.gry.gt.fnpry+iglobn-1) then print*,'Model grid point latitude must be' print*,'at least 1 gridpoint inside' print*,'pressure data area' print*,'ifm,i,j,glat,xswlat,xnelat' & ,ifm,i,j,gglat,xswlat,xnelat stop 'isan 2: outside p-data 1' endif if(idatelin.eq.1.and.gglon.lt.xswlon) gglon=gglon+360. grx=(gglon-xswlon)/gdatdx+1. if(grx.lt.2-iglobew.or.grx.gt.fnprx+iglobew-1) then print*,'Model grid point longitude must be' print*,'at least 1 gridpoint inside' print*,'pressure data area' print*,'ifm,i,j,glon,xswlon,xnelon' & ,ifm,i,j,gglon,xswlon,xnelon,iglobew,iglobs,iglobn stop 'isan 2: outside p-data 2' endif enddo enddo enddo endif ! Deallocate memory for the pressure data if(allocated(p_u)) deallocate(p_u) if(allocated(p_v)) deallocate(p_v) if(allocated(p_t)) deallocate(p_t) if(allocated(p_z)) deallocate(p_z) if(allocated(p_r)) deallocate(p_r) if(allocated(p_ur)) deallocate(p_ur) if(allocated(p_vr)) deallocate(p_vr) if(allocated(p_lat)) deallocate(p_lat) if(allocated(p_lon)) deallocate(p_lon) if(allocated(p_slp)) deallocate(p_slp) if(allocated(p_sfp)) deallocate(p_sfp) if(allocated(p_sft)) deallocate(p_sft) if(allocated(p_snow)) deallocate(p_snow) if(allocated(p_sst)) deallocate(p_sst) if(allocated(p_sc)) deallocate(p_sc) ! Allocate memory for the pressure data allocate(p_u(nprx,npry,nprz)) allocate(p_v(nprx,npry,nprz)) allocate(p_t(nprx,npry,nprz)) allocate(p_z(nprx,npry,nprz)) allocate(p_r(nprx,npry,nprz)) if(CHEM_ASSIM == 1 .and. nspecies>0 ) allocate(p_sc(nprx,npry,nprz,nspecies)) ! p_ur,p_vr arrays for rotated winds used later allocate(p_ur(nprx,npry,nprz)) allocate(p_vr(nprx,npry,nprz)) allocate(p_lat(nprx,npry)) allocate(p_lon(nprx,npry)) allocate(p_slp(nprx,npry)) allocate(p_sfp(nprx,npry)) allocate(p_sft(nprx,npry)) allocate(p_snow(nprx,npry)) allocate(p_sst(nprx,npry)) ! Fill with missing in case they are not present p_slp (1:nprx,1:npry)=1.e30 p_sfp (1:nprx,1:npry)=1.e30 p_sft (1:nprx,1:npry)=1.e30 p_snow(1:nprx,1:npry)=1.e30 p_sst (1:nprx,1:npry)=1.e30 print 192 192 format(////,1x,70('*')/ & ,' Access coarse resolution pressure data',/,1X,70('*')///) !do k=1,nprz ! pnpr(k)=levpr(k)*100. !enddo ! Call routine to fill pressure arrays from the chosen dataset. call chem_get_press_grib2() ! Find max-min theta at bottom and top levels thmax=1. thmin=1000. do j=1,npry do i=1,nprx if(p_t(i,j,nprz).lt.1.e20) & thmax=max(thmax,p_t(i,j,nprz)*(p00/pnpr(nprz))**rocp) if(p_t(i,j,1).lt.1.e20) & thmin=min(thmin,p_t(i,j,1)*(p00/pnpr(1))**rocp) enddo enddo print 300,levpr(1),thmin,levpr(nprz),thmax 300 format(//,' Minimum THETA at ',I4,' mb - ',F8.2/ & ,' Maximum THETA at ',I4,' mb - ',F8.2) ! Compute lat-lon at input pressure data points select case (inproj) case (1) do j=1,npry do i=1,nprx p_lat(i,j)=xswlat+(j-1)*gdatdy p_lon(i,j)=xswlon+(i-1)*gdatdx enddo enddo case (2) print*,'lamb-con:',xswlat,xswlon,cntlat,cntlon,gdatdx,gdatdy do j=1,npry do i=1,nprx call lc2ll(cntlat,cntlon,p_lat(i,j),p_lon(i,j) & ,float(i),float(j),xswlat,xswlon,gdatdx) ! print*,i,j,p_lat(i,j),p_lon(i,j) enddo enddo case (3) print*,'polar:',cntlat,cntlon,gdatdx,gdatdy,xswlat,xswlon xswlon_east=xswlon if(xswlon<0.) xswlon_east=xswlon+360. cntlon_east=cntlon if(cntlon<0.) cntlon_east=cntlon+360. do j=1,npry do i=1,nprx call w3fb07(float(i),float(j),xswlat,xswlon_east,gdatdx,cntlon_east & ,p_lat(i,j),p_lon(i,j)) ! print*,i,j,p_lat(i,j),p_lon(i,j) enddo enddo end select #else write(*,fmt='(A)') 'To use grib 2 the model must be compiled with grib2!' stop #endif end subroutine chem_pressure_stage_grib2 subroutine chem_get_press_grib2() use isan_coms use chem_isan_coms #ifdef GRIB2 use wgrib2api #endif use dump, only: & dumpMessage implicit none include "constants.f90" ! logical, external :: checkInside integer :: ithere(maxpr,5),isfthere(5) character*4 field,idat(5) data idat/'T','R','U','V','H'/ character(len=60) :: varName(5) real,allocatable::as(:,:) character (len=200) :: metadata character (len=300) :: grid_info character (len=99) :: invline real, allocatable :: lat(:,:),lon(:,:) real, allocatable :: aux_as(:,:) integer :: i,j,k,nv,nvar,misstot,lv,n,lv_aux integer :: irec,recordLen,iCount,jCount integer :: iAbove(5),iBellow(5) real :: vMax(5),vMin(5) integer, external :: outRealSize !Allocate array for trimmer area allocate(aux_as(nprx,npry)) iAbove=0 iBellow=0 vMax=uLimit vMin=dlimit #ifdef GRIB2 allocate(lat(nxGrib,nyGrib),lon(nxGrib,nyGrib),as(nxGrib,nyGrib)) ! Read upper air fields do lv=1,nprz_grib2 !Compose name of var to search inside grib2 varName(1)=':'//trim(wind_u_varname)//trim(slevs(lv)) varName(2)=':'//trim(wind_v_varname)//trim(slevs(lv)) varName(3)=':'//trim(temperature_varname)//trim(slevs(lv)) varName(4)=':'//trim(geo_varname)//trim(slevs(lv)) varName(5)=':'//trim(ur_varname)//trim(slevs(lv)) do nvar=1,5 !+nspecies print *,trim(innpr),trim(varName(nVar)); call flush(6) iErrNumber = grb2_inq(trim(innpr),trim(innpr)//'.inv',trim(varName(nVar)) & ,data2=as, lat=lat, lon=lon) iCount=0 jCount=1 do j=1,nyGrib do i=1,nxGrib !Teste if i,j is inside trimmer area for regional model if(checkInside(lat(i,j),lon(i,j))) then ! Mapping the position inside trimmer area iCount=iCount+1 if(iCount>nprx) then iCount=1 jCount=jCount+1 endif !Copying inside trimmer area aux_as(iCount,jCount)=as(i,j)*scale_factor(nvar) call limitVariableWithCount(varName(nVar),aux_as(iCount,jCount) & ,dlimit(nVar),uLimit(nVar),iAbove(nVar),iBellow(nVar) & ,vMax(nVar),vMin(nVar)) ! endif enddo enddo !Invert pressure leves (Grib2 -> BRAMS) and fill lv_aux=nprz_grib2-lv+1 levpr(lv_aux)=levpr_grib2(lv) pnpr(lv_aux)=levpr(lv_aux)*100 !Do not process levels above z_max_level if(lv_aux>z_max_level) cycle select case (nvar) case (1) call prfill(nprx,npry,aux_as,p_u(1,1,lv_aux)) case (2) call prfill(nprx,npry,aux_as,p_v(1,1,lv_aux)) case (3) call prfill(nprx,npry,aux_as,p_t(1,1,lv_aux)) case (4) call prfill(nprx,npry,aux_as,p_z(1,1,lv_aux)) case (5) call prfill(nprx,npry,aux_as,p_r(1,1,lv_aux)) end select print 555,int(levpr_grib2(lv)),nvar,imonth,idate,iyear,ihour 555 format(' == Read pressure field ',2I4,' at ',I2,'/',I2 & ,'/',I4,I6.4,' UTC') enddo enddo do i=1,5 !Print information about values bellow and above limits. call dumpLimits(iAbove(i),iBellow(i),varName(i),uLimit(i),dLimit(i),vMin(i),vMax(i)) enddo if(ccGradsWrite==1) call writeGradsSub() !Fill surface fields do nvar=1,5 as=-999.0 select case (nvar) case (1) call prfill(nprx,npry,as,p_slp(1,1)) case (2) call prfill(nprx,npry,as,p_sfp(1,1)) case (3) call prfill(nprx,npry,as,p_sft(1,1)) case (4) call prfill(nprx,npry,as,p_snow(1,1)) case (5) call prfill(nprx,npry,as,p_sst(1,1)) end select enddo ! Check for levels that may be all missing ithere(1:nprz,1:5)=0 isfthere(1:5)=0 do k=1,nprz do j=1,npry do i=1,nprx if(p_t(i,j,k) > 1e20) ithere(k,1)=ithere(k,1)+1 if(p_r(i,j,k) > 1e20) ithere(k,2)=ithere(k,2)+1 if(p_u(i,j,k) > 1e20) ithere(k,3)=ithere(k,3)+1 if(p_v(i,j,k) > 1e20) ithere(k,4)=ithere(k,4)+1 if(p_z(i,j,k) > 1e20) ithere(k,5)=ithere(k,5)+1 enddo enddo enddo do nv=1,5 do k=1,nprz if(ithere(k,nv) < nprx*npry) then ithere(k,nv)=1 else ithere(k,nv)=0 endif enddo enddo misstot=0 do nv=1,5 do k=1,nprz if(ithere(k,nv) == 0) misstot=misstot+1 enddo enddo do j=1,npry do i=1,nprx if(p_slp (i,j) > 1e20) isfthere(1)=isfthere(1)+1 if(p_sfp (i,j) > 1e20) isfthere(2)=isfthere(2)+1 if(p_sft (i,j) > 1e20) isfthere(3)=isfthere(3)+1 if(p_snow(i,j) > 1e20) isfthere(4)=isfthere(4)+1 if(p_sst (i,j) > 1e20) isfthere(5)=isfthere(5)+1 enddo enddo do nv=1,5 if(isfthere(nv) < nprx*npry) then isfthere(nv)=1 else isfthere(nv)=0 endif enddo print*,'------------------------------------------------' print*,' Missing parameters: 0 = all missing' print*,'------------------------------------------------' print '(t20,5(a1,6x))',(idat(n),n=1,5) do k=1,nprz print '(f10.1,t14,5(i7))',pnpr(k),(ithere(k,n),n=1,5) enddo print*,'------------------------------------------------' print*,' Missing surface parameters: 0 = all missing' print*,'------------------------------------------------' print*,' SLP SFP SFT SNOW SST' print '(t14,5(i7))',(isfthere(n),n=1,5) if(misstot.gt.0) then ! Let's see if we can get creative and make up data for the missing fields call press_miss(nprx,npry,nprz,p_u,p_v,p_t,p_z,p_r & ,ithere,maxpr,levpr) ! Check again for missing fields misstot=0 do nv=1,5 do k=1,nprz if(ithere(k,nv).eq.0) misstot=misstot+1 enddo enddo print*,'------------------------------------------------' print*,' After missing parameters check: 0 = all missing' print*,'------------------------------------------------' print '(t20,5(a1,6x))',(idat(n),n=1,5) do k=1,nprz print '(f10.1,t14,5(i7))',pnpr(k),(ithere(k,n),n=1,5) enddo endif deallocate(slevs) deallocate(mask) !write(*,fmt='(A,I4.4,2(I2.2),I4.4)') 'All done for ',iyear,imonth,idate,ihour return #endif !--(DMK-CCATT-INI)---------------------------------------------------------------- end subroutine chem_get_press_grib2 !--(DMK-CCATT-INI)---------------------------------------------------------------- subroutine chem_get_press (iunit) !--(DMK-CCATT-OLD)---------------------------------------------------------------- !subroutine get_press (iunit) !--(DMK-CCATT-FIM)---------------------------------------------------------------- use isan_coms !--(DMK-CCATT-INI)---------------------------------------------------------------- !srf-chem use chem_isan_coms !!!! use chem1_list, only : spc_name !srf-chem-end !--(DMK-CCATT-FIM)---------------------------------------------------------------- implicit none integer :: iunit real,allocatable::as(:,:) integer :: ithere(maxpr,5),isfthere(5) character*4 field,idat(5) data idat/'T','R','U','V','H'/ integer :: i,j,k,nv,nvar,misstot,lv,n integer :: recordLen,irec ! Allocate array to hold one level data for one variable allocate(as(nprx,npry)) ! Read upper air fields do lv=1,nprz !print*," ===" ,lv,nprz,nspecies;call flush(6) !--(DMK-CCATT-INI)---------------------------------------------------------------- !srf-chem !do nvar=1,5 do nvar=1,5+nspecies !srf-chem-end !srf-chem !- define the type of format that DP was writed and and reads the data if(trim(innpr(len_trim(innpr)-2:len_trim(innpr))) == 'vfm') then !print*,'reading vfm file' call vfirec(iunit,as,nprx*npry,'LIN') else !read(iunit,*,end=70,err=70) ((as(i,j),i=1,nprx),j=1,npry) print*,'reading text file *********' read(iunit,*,end=70,err=70) ((as(i,j),i=1,nprx),j=1,npry) !if(lv==1 .and. nvar==3) then ! do i=1,nprx ! do j=1,npry ! write (51,fmt='(3(F18.2,1X))') real(i),real(j) & ! ,as(i,j) ! enddo ! enddo !endif endif ! if(nvar.eq.1) then call prfill(nprx,npry,as,p_u(1,1,lv)) elseif(nvar.eq.2) then call prfill(nprx,npry,as,p_v(1,1,lv)) elseif(nvar.eq.3) then call prfill(nprx,npry,as,p_t(1,1,lv)) elseif(nvar.eq.4) then call prfill(nprx,npry,as,p_z(1,1,lv)) elseif(nvar.eq.5) then call prfill(nprx,npry,as,p_r(1,1,lv)) !--(DMK-CCATT-INI)---------------------------------------------------------------- !srf-chem elseif(nvar > 5) then !do j=1,npry;do i=1,nprx; if(as(i,j)>0.01) print*,'as=',as(i,j);enddo;enddo call prfill(nprx,npry,as,p_sc(1,1,lv,nvar-5)) !srf-chem !--(DMK-CCATT-FIM)---------------------------------------------------------------- endif print 555,levpr(lv),nvar,imonth,idate,iyear,ihour 555 format(' == Read pressure field ',2I4,' at ',I2,'/',I2 & ,'/',I4,I6.4,' UTC') enddo enddo if(ccGradsWrite==1) call writeGradsSub() ! Read surface fields do nvar=1,5 read(iunit,*,end=70,err=70) ((as(i,j),i=1,nprx),j=1,npry) ! print *,'MaxMin(as): ',maxval(as),minval(as) select case (nvar) case (1) call prfill(nprx,npry,as,p_slp(1,1)) case (2) call prfill(nprx,npry,as,p_sfp(1,1)) case (3) call prfill(nprx,npry,as,p_sft(1,1)) case (4) call prfill(nprx,npry,as,p_snow(1,1)) case (5) call prfill(nprx,npry,as,p_sst(1,1)) end select enddo goto 71 70 continue print*,'Premature end of file or error in pressure input file!' print*,'We''ll close our eyes and pretend it didn''t happen!' 71 continue deallocate(as) ! Check for levels that may be all missing ithere(1:nprz,1:5)=0 isfthere(1:5)=0 do k=1,nprz do j=1,npry do i=1,nprx if(p_t(i,j,k) > 1e20) ithere(k,1)=ithere(k,1)+1 if(p_r(i,j,k) > 1e20) ithere(k,2)=ithere(k,2)+1 if(p_u(i,j,k) > 1e20) ithere(k,3)=ithere(k,3)+1 if(p_v(i,j,k) > 1e20) ithere(k,4)=ithere(k,4)+1 if(p_z(i,j,k) > 1e20) ithere(k,5)=ithere(k,5)+1 enddo enddo enddo do nv=1,5 do k=1,nprz if(ithere(k,nv) < nprx*npry) then ithere(k,nv)=1 else ithere(k,nv)=0 endif enddo enddo misstot=0 do nv=1,5 do k=1,nprz if(ithere(k,nv) == 0) misstot=misstot+1 enddo enddo do j=1,npry do i=1,nprx if(p_slp (i,j) > 1e20) isfthere(1)=isfthere(1)+1 if(p_sfp (i,j) > 1e20) isfthere(2)=isfthere(2)+1 if(p_sft (i,j) > 1e20) isfthere(3)=isfthere(3)+1 if(p_snow(i,j) > 1e20) isfthere(4)=isfthere(4)+1 if(p_sst (i,j) > 1e20) isfthere(5)=isfthere(5)+1 enddo enddo do nv=1,5 if(isfthere(nv) < nprx*npry) then isfthere(nv)=1 else isfthere(nv)=0 endif enddo print*,'------------------------------------------------' print*,' Missing parameters: 0 = all missing' print*,'------------------------------------------------' print '(t20,5(a1,6x))',(idat(n),n=1,5) do k=1,nprz print '(f10.1,t14,5(i7))',pnpr(k),(ithere(k,n),n=1,5) enddo print*,'------------------------------------------------' print*,' Missing surface parameters: 0 = all missing' print*,'------------------------------------------------' print*,' SLP SFP SFT SNOW SST' print '(t14,5(i7))',(isfthere(n),n=1,5) if(misstot.gt.0) then ! Let's see if we can get creative and make up data for the missing fields call press_miss(nprx,npry,nprz,p_u,p_v,p_t,p_z,p_r & ,ithere,maxpr,levpr) ! Check again for missing fields misstot=0 do nv=1,5 do k=1,nprz if(ithere(k,nv).eq.0) misstot=misstot+1 enddo enddo print*,'------------------------------------------------' print*,' After missing parameters check: 0 = all missing' print*,'------------------------------------------------' print '(t20,5(a1,6x))',(idat(n),n=1,5) do k=1,nprz print '(f10.1,t14,5(i7))',pnpr(k),(ithere(k,n),n=1,5) enddo endif return !--(DMK-CCATT-INI)---------------------------------------------------------------- end subroutine chem_get_press #ifdef netcdf subroutine chem_pressure_stage_netCDF(n1,n2,nhem,glat,glon,glat2,glon2) !--(DMK-CCATT-OLD)---------------------------------------------------------------- !subroutine pressure_stage(n1,n2,nhem,glat,glon,glat2,glon2) !--(DMK-CCATT-FIM)---------------------------------------------------------------- use isan_coms use ModDateUtils, only: & date_add_to, date_add_to_dble !--(DMK-CCATT-INI)---------------------------------------------------------------- !srf-chem use chem_isan_coms use chem1_list, only: chemical_mechanism ! intent(in) use mem_chem1, only: CHEM_ASSIM, CHEMISTRY ! intent(in) !srf-chem-end !--(DMK-CCATT-FIM)---------------------------------------------------------------- use rconstants use netcdf use dump, only: & dumpMessage implicit none include "constants.f90" include "netcdf.inc" character(len=*),parameter :: header='**(chem_pressure_stage_netCDF)**' logical, external :: checkInside integer :: n1,n2,nhem real :: glat(n1,n2),glon(n1,n2),glat2(n1,n2),glon2(n1,n2) real :: fnprx,fnpry,grx,gry,gglat,gglon,thmax,thmin & ,xswlon_east,cntlon_east,rr,each_levp_real(maxpr) integer :: i,j,k,lv,ifm,loop,n,iunit,stat real, external :: rs character (len=200) :: metadata character(len=32),allocatable :: varName(:) character(len=32) :: name,atName !# composed name of levels (:press:) character (len=300) :: grid_info character (len=99) :: invline real, allocatable :: var(:,:),lat(:,:),lon(:,:) real, allocatable :: llat(:),llon(:),pLev(:) integer :: nx,ny,thisHour,iFound,iFactor integer :: iyyl,imml,iddl,ihhl,an,iLen,xtype integer :: ndims, nvars, nglobalatts, unlimdimid,ncid integer, allocatable :: nat(:),varDim(:),lenDim(:) character(len=32),allocatable :: atValue(:) integer :: levelVarN,lonVarN,latVarN,uVarN,vVarN,tVarN,zVarN,qVarN integer :: nTimes,timeVarN,validTime integer, allocatable :: HoursFrom1900(:) real, allocatable :: latTest(:) real :: xvalue real :: scaleFactor(32),addOfset(32) integer :: ii,jj,kk integer :: iAbove(5),iBellow(5) real :: vMax(5),vMin(5) integer :: beginDate,beginTime character(len=33) :: timeLong real :: missingValue !--(DMK-CCATT-INI)---------------------------------------------------------------- !srf-chem character(len=32) :: chemical_mechanism_test !srf-chem-end !--(DMK-CCATT-FIM)---------------------------------------------------------------- call str2int(innpr(len_trim(innpr)-2:len_trim(innpr)),thisHour,stat) !Open NetCDF file !write(*,fmt='(A,A,A)') '{',trim(innpr),'}' iErrNumber = nf90_open(path = trim(innpr), mode = nf90_nowrite, ncid = ncid) if (iErrNumber /= nf90_noerr) iErrNumber=dumpMessage(c_tty,c_yes,header,modelVersion & ,c_fatal,trim(innpr)//' not found. Please, check it!') print 91,innpr(1:len_trim(innpr)) 91 format(//,' Reading pressure gridded netCDF data',/,a,//) ! get info about netCDF file iErrNumber=nf90_inquire(ncid, ndims, nvars) !print *,'File open. Total of variables are ',nvars allocate(varName(nvars),varDim(nvars),nat(nvars),lenDim(nvars)) ! Get dimensions in levels, lons, lats and times nprz_grib2=0 nxGrib=0 nyGrib=0 do i=1,nvars iErrNumber = nf90_Inquire_Dimension(ncid, i, name, lenDim(i)) !print *,'Variable ',i,'name ',trim(name) select case (trim(name)) case ('level') nprz_grib2=lenDim(i) levelVarN=i case ('lev') nprz_grib2=lenDim(i) levelVarN=i case ('longitude') nxGrib=lenDIm(i) lonVarN=i case ('lon') nxGrib=lenDIm(i) lonVarN=i case ('latitude') nyGrib=lenDIm(i) latVarN=i case ('lat') nyGrib=lenDIm(i) latVarN=i case ('time') nTimes=lenDIm(i) timeVarN=i end select enddo !If one of variables not found is error if(nprz_grib2==0 .or. nxGrib==0 .or. nyGrib==0) then print *,'Vars expected: level, longitude, latitude and time.' iErrNumber=dumpMessage(c_tty,c_yes,header,modelVersion & ,c_fatal,'Var expected for dimension not found in netCDF file. Please check it!') endif write(*,fmt='("number of lats points: ",I4)') nygrib write(*,fmt='("number of lons points: ",I4)') nxgrib write(*,fmt='("number of levs : ",I4)') nprz_grib2 !Allocate all vars allocate(PressLevs(nprz_grib2)) allocate(mask(nxGrib,nyGrib)) allocate(llat(nyGrib),llon(nxGrib)) allocate(HoursFrom1900(nTimes)) do i=1,5 allocate(ncVar(i)%valuesRead(nxGrib,nyGrib,nprz_grib2)) enddo !Get pressure levels iErrNumber=nf90_get_var(ncid, levelVarN, PressLevs) do i=1,nprz_grib2 levpr_grib2(i)=PressLevs(i) enddo !If the total of levels ir greater than defined in RAMSIN use the defined if(nprz_grib2>z_max_level) then nprz=z_max_level else nprz=nprz_grib2 endif !Get names of vars in netCDF file and total of atributes do i=1,nvars iErrNumber=nf90_Inquire_Variable(ncid, i, name=varName(i) & ,ndims=varDim(i), nAtts=nat(i)) !print *,i,trim(varName(i)),varDim(i),nat(i) enddo !checking if required varName is present iFound=0 do i=1,nvars select case (trim(varName(i))) case (trim(WIND_U_VARNAME)) iFound=iFound+1 uVarN=i case (trim(WIND_V_VARNAME)) iFound=iFound+1 vVarN=i case (trim(TEMPERATURE_VARNAME)) iFound=iFound+1 tVarN=i case (trim(GEO_VARNAME)) iFound=iFound+1 zVarN=i case (trim(UR_VARNAME)) iFound=iFound+1 qVarN=i end select enddo if(iFound/=5) then !Not found all needs print *,'Available names of var in netCDF file:' do i=1,nvars print *,trim(varName(i)) enddo print *,'Expected name os vars from RAMSIN:' print *,trim(WIND_U_VARNAME) print *,trim(WIND_V_VARNAME) print *,trim(TEMPERATURE_VARNAME) print *,trim(GEO_VARNAME) print *,trim(UR_VARNAME) iErrNumber=dumpMessage(c_tty,c_yes,header,modelVersion & ,c_fatal,'Incompatible var names. See the list above and RAMSIN vars!') endif missingValue=1.0e+20 !Getting variables atributes do i=1,nvars do an=1,nat(i) iErrNumber=nf90_inq_attname(ncid, i, an, atName) iErrNumber=nf90_Inquire_Attribute(ncid, i, trim(atName), xtype=xtype, len=iLen) if(trim(atName)=='scale_factor') then iErrNumber = nf90_get_att(ncid, i, atName, xValue) scaleFactor(i)=xvalue elseif(trim(atName)=='add_offset') then iErrNumber = nf90_get_att(ncid, i, atName, xValue) addOfset(i)=xvalue elseif(trim(atName)=='begin_date') then iErrNumber = nf90_get_att(ncid, i, atName, xValue) beginDate=xValue elseif(trim(atName)=='begin_time') then iErrNumber = nf90_get_att(ncid, i, atName, xValue) beginTime=xValue elseif(trim(atName)=='units' .and. i==timeVarN) then iErrNumber = nf90_get_att(ncid, i, atName, timeLong) elseif(trim(atName)=='missing_value') then iErrNumber = nf90_get_att(ncid, i, atName, missingValue) endif enddo end do !get lon and lat info iErrNumber = nf90_get_var(ncid, latVarN, llat) iErrNumber = nf90_get_var(ncid, lonVarN, llon) if(ICFILETYPE==3) then !NASA data is from -90 to +90 gdatdy=llat(2)-llat(1) else gdatdy=llat(1)-llat(2) !ECMWF data is from +90 to -90 endif gdatdx=llon(2)-llon(1) write(*,fmt='("First lat point : ",F8.2)') llat(1) write(*,fmt='("First lon point : ",F8.2)') llon(1) write(*,fmt='("Delta lat : ",F8.4)') gdatdy write(*,fmt='("Delta lon : ",F8.4)') gdatdx write(*,fmt='("Missing Value : ",E18.4)') missingValue nprx=((final_longitude-initial_longitude)/gdatdx)+1 npry=((final_latitude-initial_latitude)/gdatdy)+1 !Setting initial lon and lat xswlon=initial_longitude-360.0 xswlat=initial_latitude ! Check for date consistency between file parameters and namelist parameters iErrNumber = nf90_get_var(ncid, timeVarN, HoursFrom1900) if(iErrNumber /= nf90_noerr) print *, 'Err: ',iErrNumber !iErrNumber = nf90_get_var(ncid, latVarN, latTest) !if(iErrNumber /= nf90_noerr) print *, 'Err: ',iErrNumber !print *, "LatTest: ",size(latTest),latTest write (*,fmt='(A,I2,A,I2,A)') "Var time from netCDF is #",timeVarN," with ",size(HoursFrom1900),' times!' print *,HoursFrom1900,nTimes validTime=0 if(ICFILETYPE==3) then call splitDate(timeLong,iyy,imm,idd,ihh) !print *,timelong,iyy,imm,idd,ihh print *,'Time: ',iyy,imm,idd,ihh,iyear,imonth,idate,ihour if(iyy==iyear .and. imm==imonth .and. idd==idate & .and. ihh==ihour) then validTime=1 endif else do i=1,nTimes write (*,fmt='(A,I2,A,I16)') "Hours for ",i," is ",HoursFrom1900(i) call date_add_to_dble(1900,1,1,0,dble(HoursFrom1900(i)),'h' & ,iyy,imm,idd,ihh) if(iyy==iyear .and. imm==imonth .and. idd==idate & .and. ihh/100==ihour) then validTime=i exit endif enddo endif if(validTime==0) then !Not found the date inside netCDF print*,'Pressure file dates not the same as namelist!' print *,'Namelist require: ',iyear,imonth,idate,ihour print *,'The netCDF file has:' do i=1,nTimes call date_add_to_dble(1900,1,1,0,dble(HoursFrom1900(i)),'h' & ,iyy,imm,idd,ihh) write (*,fmt='(I2.2,1X,I4.4,1X,I2.2,1X,I2.2,1X,I8)') i,iyy,imm,idd,ihh,HoursFrom1900(i) enddo iErrNumber=dumpMessage(c_tty,c_yes,header,modelVersion & ,c_fatal,'Please, check the dates!') endif !Read all upper air values for u,v,t,z and r do i=1,5 select case (i) case(1) iErrNumber = nf90_get_var(ncid, uVarN, ncVar(i)%valuesRead & ,start = (/ 1, 1, 1,validTime /)) call adjustMissingValues(icFileType,i,nxGrib,nyGrib,nprz_grib2,missingValue & ,levpr_grib2,ncVar(i)%valuesRead) ncVar(i)%valuesRead=ncVar(i)%valuesRead*scaleFactor(uVarN) & +addOfset(uVarN) case(2) iErrNumber = nf90_get_var(ncid, vVarN, ncVar(i)%valuesRead & ,start = (/ 1, 1, 1,validTime /)) call adjustMissingValues(icFileType,i,nxGrib,nyGrib,nprz_grib2,missingValue & ,levpr_grib2,ncVar(i)%valuesRead) ncVar(i)%valuesRead=ncVar(i)%valuesRead*scaleFactor(vVarN) & +addOfset(vVarN) case(3) iErrNumber = nf90_get_var(ncid, tVarN, ncVar(i)%valuesRead & ,start = (/ 1, 1, 1,validTime /)) call adjustMissingValues(icFileType,i,nxGrib,nyGrib,nprz_grib2,missingValue & ,levpr_grib2,ncVar(i)%valuesRead) ncVar(i)%valuesRead=ncVar(i)%valuesRead*scaleFactor(tVarN) & +addOfset(tVarN) case(4) iErrNumber = nf90_get_var(ncid, zVarN, ncVar(i)%valuesRead & ,start = (/ 1, 1, 1,validTime /)) call adjustMissingValues(icFileType,i,nxGrib,nyGrib,nprz_grib2,missingValue & ,levpr_grib2,ncVar(i)%valuesRead) ncVar(i)%valuesRead=ncVar(i)%valuesRead*scaleFactor(zVarN) & +addOfset(zVarN) case(5) iErrNumber = nf90_get_var(ncid, qVarN, ncVar(i)%valuesRead & ,start = (/ 1, 1, 1,validTime /)) call adjustMissingValues(icFileType,i,nxGrib,nyGrib,nprz_grib2,missingValue & ,levpr_grib2,ncVar(i)%valuesRead) ncVar(i)%valuesRead=ncVar(i)%valuesRead*scaleFactor(qVarN) & +addOfset(qVarN) end select !Adjust the value accordingly RAMSIN scale_factor ncVar(i)%valuesRead=ncVar(i)%valuesRead*scale_factor(i) !If values out of valid values setup to limits iAbove(i)=0 iBellow(i)=0 vMax(i)=uLimit(i) vMin(i)=dlimit(i) do ii=1,nxGrib do jj=1,nyGrib if(.not. checkInside(llat(jj),llon(ii))) cycle do kk=1,nprz_grib2 call limitVariableWithCount(varName(i),ncVar(i)%valuesRead(ii,jj,kk) & ,dlimit(i),uLimit(i),iAbove(i),iBellow(i),vMax(i),vMin(i)) enddo enddo enddo !Print information about values bellow and above limits. call dumpLimits(iAbove(i),iBellow(i),varName(i),uLimit(i),dLimit(i),vMin(i),vMax(i)) enddo if(ICFILETYPE==3) then inproj=1 !NASA lons is from -180 to +180 else inproj=1 !ECMWF lons is from 0 to 360 endif ! Check pressure data domain size and location if (inproj/=1.and.nhem>0) then print*,'You must input a lat-lon pressure grid ' & ,'to run a global simulation !!!!' stop 'glob-no-press' endif if (inproj.eq.1) then ! If necessary, convert longitude specification to values in the range ! [-180.,180.]. xswlon=mod(xswlon+900.,360.)-180. xnelon=mod(xnelon-900.,-360.)+180. if(xswlon.lt.-180..or.xswlon.ge.180. .or. & xnelon.lt.-180..or.xnelon.gt.180.01.or. & xswlat.lt.-90. .or.xnelat.gt.90.01) then print*,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' print*,'!!! BAD DOMAIN SPECIFICATION !!' print*,'!!! xswlat,xswlon - ',xswlat,xswlon print*,'!!! xnelat,xnelon-xswlon - ',xnelat,xnelon-xswlon print*,'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' stop 'bad_domain' endif fnprx=float(nprx) fnpry=float(npry) ! Set "global domain" flags to determine whether 2 extra rows will be ! added to borders of the input gridded data. iglobew = 1 if doing ! full 360 degrees of longitude. iglobs = 1 if including the south ! pole and if iglobew = 1, and iglobn = 1 if including ! the north pole and if iglobew = 1. iglobew=0 iglobn=0 iglobs=0 if(fnprx*gdatdx.gt.359.9) then if(abs(xswlon+180.).gt.1.e-3) then print*,'When extracting all longitudes of pressure data,' print*,'xswlon must be set to -180.' !stop 'xswlon' endif iglobew=1 endif if(iglobew.eq.1.and.xswlat.lt.-89.9) iglobs=1 if(iglobew.eq.1.and.xnelat.gt. 89.9) iglobn=1 idatelin=0 if(xnelon.lt.xswlon+.01) idatelin=1 ! Search coarse grid points for any that are outside the bounds of ! the pressure data. If any are found, stop. loop=1 if(nhem>0) loop=2 do ifm = 1,loop do j = 1,n2 do i = 1,n1 if(ifm==1) then gglat=glat(i,j) gglon=glon(i,j) elseif(ifm==2) then gglat=glat2(i,j) gglon=glon2(i,j) endif gry = (gglat - xswlat) / gdatdy + 1. if(gry .lt. 2-iglobs .or. gry .gt. fnpry+iglobn-1) then write(*,fmt='("i : ",I4.4)') i write(*,fmt='("j : ",I4.4)') j write(*,fmt='("xswlat : ",E18.4)') xswlat write(*,fmt='("xnelat : ",E18.4)') xnelat write(*,fmt='("gLat : ",E18.4)') gglat write(*,fmt='("xnelat : ",E18.4)') xnelat write(*,fmt='("iglobs : ", I1)') iglobs write(*,fmt='("iglobn : ", I1)') iglobn write(*,fmt='("fnpry : ",E18.4)') fnpry write(*,fmt='("gdatdy : ",E18.4)') gdatdy write(*,fmt='("gry : ",E18.4)') gry if(gry.lt.2-iglobs) write(*,fmt='(A)') "gry is < 2-iglobs!" if(gry.gt.fnpry+iglobn-1) write(*,fmt='(A)') "gry is > fnpry+iglobn-1!" iErrNumber=dumpMessage(c_tty,c_yes,header,'Pos 3' & ,c_fatal,'Model grid point latitude must be at least 1' & //'gridpoint inside pressure data area, see info above!') endif if(idatelin.eq.1.and.gglon.lt.xswlon) gglon=gglon+360. grx=(gglon-xswlon)/gdatdx+1. if(grx.lt.2-iglobew.or.grx.gt.fnprx+iglobew-1) then print*,'Model grid point longitude must be' print*,'at least 1 gridpoint inside' print*,'pressure data area' print*,'ifm,i,j,glon,xswlon,xnelon' & ,ifm,i,j,gglon,xswlon,xnelon,iglobew,iglobs,iglobn stop 'isan 3: outside p-data 2' endif enddo enddo enddo endif ! Deallocate memory for the pressure data if(allocated(p_u)) deallocate(p_u) if(allocated(p_v)) deallocate(p_v) if(allocated(p_t)) deallocate(p_t) if(allocated(p_z)) deallocate(p_z) if(allocated(p_r)) deallocate(p_r) if(allocated(p_ur)) deallocate(p_ur) if(allocated(p_vr)) deallocate(p_vr) if(allocated(p_lat)) deallocate(p_lat) if(allocated(p_lon)) deallocate(p_lon) if(allocated(p_slp)) deallocate(p_slp) if(allocated(p_sfp)) deallocate(p_sfp) if(allocated(p_sft)) deallocate(p_sft) if(allocated(p_snow)) deallocate(p_snow) if(allocated(p_sst)) deallocate(p_sst) if(allocated(p_sc)) deallocate(p_sc) ! Allocate memory for the pressure data allocate(p_u(nprx,npry,nprz)) allocate(p_v(nprx,npry,nprz)) allocate(p_t(nprx,npry,nprz)) allocate(p_z(nprx,npry,nprz)) allocate(p_r(nprx,npry,nprz)) if(CHEM_ASSIM == 1 .and. nspecies>0 ) allocate(p_sc(nprx,npry,nprz,nspecies)) ! p_ur,p_vr arrays for rotated winds used later allocate(p_ur(nprx,npry,nprz)) allocate(p_vr(nprx,npry,nprz)) allocate(p_lat(nprx,npry)) allocate(p_lon(nprx,npry)) allocate(p_slp(nprx,npry)) allocate(p_sfp(nprx,npry)) allocate(p_sft(nprx,npry)) allocate(p_snow(nprx,npry)) allocate(p_sst(nprx,npry)) ! Fill with missing in case they are not present p_slp (1:nprx,1:npry)=1.e30 p_sfp (1:nprx,1:npry)=1.e30 p_sft (1:nprx,1:npry)=1.e30 p_snow(1:nprx,1:npry)=1.e30 p_sst (1:nprx,1:npry)=1.e30 print 192 192 format(////,1x,70('*')/ & ,' Access coarse resolution pressure data',/,1X,70('*')///) !Call routine to fill pressure arrays from the chosen dataset. call chem_get_press_netCDF(uVarN,vVarN,tVarN,zVarN,qVarN & ,llat,llon,nTimes) ! Find max-min theta at bottom and top levels thmax=1. thmin=1000. do j=1,npry do i=1,nprx if(p_t(i,j,nprz).lt.1.e20) & thmax=max(thmax,p_t(i,j,nprz)*(p00/pnpr(nprz))**rocp) if(p_t(i,j,1).lt.1.e20) & thmin=min(thmin,p_t(i,j,1)*(p00/pnpr(1))**rocp) enddo enddo print 300,levpr(1),thmin,levpr(nprz),thmax 300 format(//,' Minimum THETA at ',I4,' mb - ',F8.2/ & ,' Maximum THETA at ',I4,' mb - ',F8.2) ! Compute lat-lon at input pressure data points select case (inproj) case (1) do j=1,npry do i=1,nprx p_lat(i,j)=xswlat+(j-1)*gdatdy p_lon(i,j)=xswlon+(i-1)*gdatdx enddo enddo case (2) print*,'lamb-con:',xswlat,xswlon,cntlat,cntlon,gdatdx,gdatdy do j=1,npry do i=1,nprx call lc2ll(cntlat,cntlon,p_lat(i,j),p_lon(i,j) & ,float(i),float(j),xswlat,xswlon,gdatdx) ! print*,i,j,p_lat(i,j),p_lon(i,j) enddo enddo case (3) print*,'polar:',cntlat,cntlon,gdatdx,gdatdy,xswlat,xswlon xswlon_east=xswlon if(xswlon<0.) xswlon_east=xswlon+360. cntlon_east=cntlon if(cntlon<0.) cntlon_east=cntlon+360. do j=1,npry do i=1,nprx call w3fb07(float(i),float(j),xswlat,xswlon_east,gdatdx,cntlon_east & ,p_lat(i,j),p_lon(i,j)) ! print*,i,j,p_lat(i,j),p_lon(i,j) enddo enddo end select iErrNumber = nf90_close(ncid) end subroutine chem_pressure_stage_netCDF subroutine chem_get_press_netCDF(uVarN,vVarN,tVarN,zVarN,qVarN & ,llat,llon,nTimes) use isan_coms use chem_isan_coms use dump, only: & dumpMessage use netcdf implicit none include "constants.f90" ! character(len=*), parameter :: cFMT='(" == Read pressure field ",I4," mBar for var=",A1," at ",I2,"/",I2,"/",I4,I6.4," UTC, maxval= ",E18.6,", minval= ",E18.6)' integer,intent(in) :: uVarN,vVarN,tVarN,zVarN,qVarN,nTimes real,intent(in) :: llat(nyGrib),llon(nxGrib) integer :: ncid logical, external :: checkInside integer :: ithere(maxpr,5),isfthere(5) character*4 field,idat(5) data idat/'T','R','U','V','H'/ character,parameter,dimension(5) :: varName=(/'U','V','T','Z','R'/) real,allocatable::as(:,:) character (len=200) :: metadata character (len=300) :: grid_info character (len=99) :: invline real, allocatable :: lat(:,:),lon(:,:) real, allocatable :: aux_as(:,:),aux_as_(:,:) integer :: i,j,k,nv,nvar,misstot,lv,n,lv_aux,auxJ integer :: irec,recordLen,iCount,jCount,iFactor !Allocate array for trimmer area allocate(aux_as(nprx,npry)) allocate(aux_as_(nprx,npry)) allocate(as(nxGrib,nyGrib)) ! Put the upper air info using correct lats and lons and invert levels do lv=nprz_grib2,1,-1 do nvar=1,5 as=ncVar(nvar)%valuesRead(:,:,lv) !if(nvar==3) write(55,fmt='(I2.2,1X,6(F8.2,1X),10(E18.3,1X))') lv & ! ,llat(151),llon(216) & ! ,initial_latitude,final_latitude,initial_longitude,final_longitude & ! ,as(251,161),as(216,151) & ! ,as(251,162),as(216,152) & ! ,as(251,163),as(216,153) & ! ,as(251,164),as(216,154) & ! ,as(251,165),as(216,155) iCount=0 jCount=1 do j=nyGrib,1,-1 do i=1,nxGrib !Teste if i,j is inside trimmer area for regional model if(checkInside(llat(j),llon(i))) then ! Mapping the position inside trimmer area iCount=iCount+1 if(iCount>nprx) then iCount=1 jCount=jCount+1 endif !Copying inside trimmer area aux_as(iCount,jCount)=as(i,j) ! if(nvar==3 .and. i==216 .and. j==151) write(55,fmt='(2(I2.2,1X),E18.3)') & ! icount,jcount,aux_as(iCount,jCount) endif enddo enddo if(ICFILETYPE==3) then !Inverting Latitudes (NASA is from -90 to +90) lv_aux=lv do j=jCount,1,-1 auxJ=jCount-j+1 aux_as_(:,j)=aux_as(:,auxJ) enddo aux_as=aux_as_ else !Inverting pressure levels (netCDF -> BRAMS) and fill lv_aux=nprz_grib2-lv+1 !Others is from up to down endif levpr(lv_aux)=levpr_grib2(lv) pnpr(lv_aux)=levpr(lv_aux)*100 ! !Do not process levels above z_max_level if(lv_aux>z_max_level) cycle select case (nvar) case (1) call prfill(nprx,npry,aux_as,p_u(1,1,lv_aux)) case (2) call prfill(nprx,npry,aux_as,p_v(1,1,lv_aux)) case (3) call prfill(nprx,npry,aux_as,p_t(1,1,lv_aux)) case (4) call prfill(nprx,npry,aux_as,p_z(1,1,lv_aux)) case (5) call prfill(nprx,npry,aux_as,p_r(1,1,lv_aux)) end select write(*,fmt=trim(cFMT)) int(levpr_grib2(lv)),varName(nvar),imonth,idate,iyear,ihour & ,maxval(aux_as),minval(aux_as) enddo enddo if(ccGradsWrite==1) call writeGradsSub() !Fill surface fields do nvar=1,5 as=-999.0 select case (nvar) case (1) call prfill(nprx,npry,as,p_slp(1,1)) case (2) call prfill(nprx,npry,as,p_sfp(1,1)) case (3) call prfill(nprx,npry,as,p_sft(1,1)) case (4) call prfill(nprx,npry,as,p_snow(1,1)) case (5) call prfill(nprx,npry,as,p_sst(1,1)) end select enddo ! Check for levels that may be all missing ithere(1:nprz,1:5)=0 isfthere(1:5)=0 do k=1,nprz do j=1,npry do i=1,nprx if(p_t(i,j,k) > 1e20) ithere(k,1)=ithere(k,1)+1 if(p_r(i,j,k) > 1e20) ithere(k,2)=ithere(k,2)+1 if(p_u(i,j,k) > 1e20) ithere(k,3)=ithere(k,3)+1 if(p_v(i,j,k) > 1e20) ithere(k,4)=ithere(k,4)+1 if(p_z(i,j,k) > 1e20) ithere(k,5)=ithere(k,5)+1 enddo enddo enddo do nv=1,5 do k=1,nprz if(ithere(k,nv) < nprx*npry) then ithere(k,nv)=1 else ithere(k,nv)=0 endif enddo enddo misstot=0 do nv=1,5 do k=1,nprz if(ithere(k,nv) == 0) misstot=misstot+1 enddo enddo do j=1,npry do i=1,nprx if(p_slp (i,j) > 1e20) isfthere(1)=isfthere(1)+1 if(p_sfp (i,j) > 1e20) isfthere(2)=isfthere(2)+1 if(p_sft (i,j) > 1e20) isfthere(3)=isfthere(3)+1 if(p_snow(i,j) > 1e20) isfthere(4)=isfthere(4)+1 if(p_sst (i,j) > 1e20) isfthere(5)=isfthere(5)+1 enddo enddo do nv=1,5 if(isfthere(nv) < nprx*npry) then isfthere(nv)=1 else isfthere(nv)=0 endif enddo print*,'------------------------------------------------' print*,' Missing parameters: 0 = all missing' print*,'------------------------------------------------' print '(t20,5(a1,6x))',(idat(n),n=1,5) do k=1,nprz print '(f10.1,t14,5(i7))',pnpr(k),(ithere(k,n),n=1,5) enddo print*,'------------------------------------------------' print*,' Missing surface parameters: 0 = all missing' print*,'------------------------------------------------' print*,' SLP SFP SFT SNOW SST' print '(t14,5(i7))',(isfthere(n),n=1,5) if(misstot.gt.0) then ! Let's see if we can get creative and make up data for the missing fields call press_miss(nprx,npry,nprz,p_u,p_v,p_t,p_z,p_r & ,ithere,maxpr,levpr) ! Check again for missing fields misstot=0 do nv=1,5 do k=1,nprz if(ithere(k,nv).eq.0) misstot=misstot+1 enddo enddo print*,'------------------------------------------------' print*,' After missing parameters check: 0 = all missing' print*,'------------------------------------------------' print '(t20,5(a1,6x))',(idat(n),n=1,5) do k=1,nprz print '(f10.1,t14,5(i7))',pnpr(k),(ithere(k,n),n=1,5) enddo endif deallocate(mask) deallocate(pressLevs) do nVar=1,5 deallocate(ncVar(nvar)%valuesRead) enddo !write(*,fmt='(A,I4.4,2(I2.2),I4.4)') 'All done for ',iyear,imonth,idate,ihour return end subroutine chem_get_press_netCDF #endif !!=============================== !Utils functions and subroutines logical function checkInside(xlat,xlon) use isan_coms real, intent(in) :: xlat,xlon checkInside=.false. if(xlat>=initial_latitude .and. xlat<=final_latitude & .and. xlon>=initial_longitude .and. xlon & <=final_longitude) checkInside=.true. end function checkInside subroutine limitVariableWithCount(variableName, variable, inferiorLim & , superiorLim,above,bellow,vMax,vMin) use dump, only: & dumpMessage include "constants.f90" character(len=*),parameter :: header='**(limitVariable)**' real, intent(inout) :: variable integer, intent(inout) :: above integer, intent(inout) :: bellow real, intent(inout) :: vMax real, intent(inout) :: vMin real, intent(in) :: inferiorLim, superiorLim character(len=*), intent(in) :: variableName if (variable .lt. inferiorLim) then if(variablevMax) vmax=variable variable=superiorLim above=above+1 endif end subroutine limitVariableWithCount subroutine dumpLimits(iAbove,iBellow,varName,uLimit,dLimit,vMin,vMax) !Header_vai_aqui! use dump, only: & dumpMessage implicit none include "constants.f90" !Parameters (constants) ! Input/Output variables integer,intent(in) :: iAbove,iBellow !# real ,intent(in) :: uLimit,dLimit,vMin,vMax !# character(len=*),intent(in) :: varName !# !Code if(iAbove>0) then iErrNumber=dumpMessage(c_tty,c_yes,'','',c_warning & ,'The maximum value allowed for '//trim(varName) & //' is:',uLimit,'E18.4') iErrNumber=dumpMessage(c_tty,c_yes,'','',c_warning & ,'There are values above the limit ' & //' - total count:',iabove,'I8') iErrNumber=dumpMessage(c_tty,c_yes,'','',c_warning & ,'The greater value found is:',vMax,'E18.4') endif if(iBellow>0) then iErrNumber=dumpMessage(c_tty,c_yes,'','',c_warning & ,'The minimum value allowed for '//trim(varName) & //' is:',dLimit,'E18.4') iErrNumber=dumpMessage(c_tty,c_yes,'','',c_warning & ,'There are values bellow the limit ' & //' - total count:',iBellow,'I8') iErrNumber=dumpMessage(c_tty,c_yes,'','',c_warning & ,'The minor value is:',vMin,'E18.4') endif end subroutine dumpLimits subroutine writeGradsSub() !Header_vai_aqui! use isan_coms implicit none include "constants.f90" !Parameters (constants) ! Input/Output variables !Local variables integer :: recordLen, irec, nvar,k character(len=4) :: cYear,cHour character(len=2) :: cDay,cMonth character(len=256) :: gradsFileName character(len=19) :: fileName integer, external :: outRealSize !Code write(cYear,fmt='(I4.4)') iyear write(cMonth,fmt='(I2.2)') iMonth write(cDay,fmt='(I2.2)') idate write(cHour,fmt='(I4.4)') ihour fileName='ciREad.'//cYear//cMonth//cDay//cHour gradsFileName=trim(ICGRADSPREFIX)//fileName recordLen=outRealSize()*nprx*npry ! recordLen=nprx*npry open(unit=33,file=trim(gradsFileName)//'.gra',& action='WRITE',status='REPLACE',form='UNFORMATTED',access='DIRECT', & recl=recordLen) !# writing grads binary and fill variables irec=1 do nvar=1,5 do k=1,nprz select case (nvar) case (1) write (33,rec=irec) p_u (1:nprx,1:npry,k) case (2) write (33,rec=irec) p_v (1:nprx,1:npry,k) case (3) write (33,rec=irec) p_t (1:nprx,1:npry,k) case (4) write (33,rec=irec) p_z (1:nprx,1:npry,k) case (5) write (33,rec=irec) p_r (1:nprx,1:npry,k) end select irec=irec+1 enddo enddo close(33) open(unit=33,file=trim(gradsFileName)//'.ctl' & ,action='WRITE',status='replace',form='FORMATTED') !writing the name of grads file write(33,*) 'dset ^'//fileName//'.gra' !writing others infos to ctl write(33,*) 'undef -0.9990000E+34' write(33,*) 'title GRIB NCEP GFS' write(33,*) 'xdef ',nprx,' linear ',xswlon,gdatdx write(33,*) 'ydef ',npry,' linear ',xswlat,gdatdy write(33,*) 'zdef ',nprz,'levels',levpr(1:nprz) write(33,*) 'tdef 1 linear 00:00z01jan2018 1mo' write(33,*) 'vars 5' write(33,*) 'U ',nprz,'99 U wind' write(33,*) 'V ',nprz,'99 V Wind' write(33,*) 'T ',nprz,'99 Temp' write(33,*) 'Z ',nprz,'99 Geopot' write(33,*) 'R ',nprz,'99 Rel Umid' write(33,*) 'endvars' close(33) end subroutine writeGradsSub subroutine str2int(str,int,stat) implicit none ! Arguments character(len=*),intent(in) :: str integer,intent(out) :: int integer,intent(out) :: stat read(str,*,iostat=stat) int end subroutine str2int subroutine limitVariable(variableName, variable, inferiorLim, superiorLim) use dump, only: & dumpMessage include "constants.f90" character(len=*),parameter :: header='**(limitVariable)**' real, intent(inout) :: variable real, intent(in) :: inferiorLim, superiorLim character(len=*), intent(in) :: variableName if (variable .lt. inferiorLim) then write(12,fmt='("Variable, Limit: ",2(F18.2,1X))') variable,inferiorLim iErrNumber=dumpMessage(12,c_yes,header,modelVersion & ,c_warning,'Value of '& //trim(variableName)//' is inferior than limit.' & //' Using limit.') variable=inferiorLim else if ( variable .gt. superiorLim) then write(12,fmt='("Variable, Limit: ",2(F18.2,1X))') variable,superiorLim iErrNumber=dumpMessage(12,c_yes,header,modelVersion & ,c_warning,'Value of '& //trim(variableName)//' is superior than limit.' & //' Using limit.') variable=superiorLim endif end subroutine limitVariable subroutine splitDate(di,yy,mm,dd,hh) !# Subroutine to convert GEOS NetCDF date long unit to !# year, month, day and hour !# !# @note !# ![](http://brams.cptec.inpe.br/wp-content/uploads/2015/11/logo-brams.navigation.png "") !# !# **Brief**: This subroutine trim the lond unit date description !# minutes since 2020-02-13 00:00:00 into year, month, day and hour !# !# **Documentation**: !# !# **Author**: Luiz Flavio Rodrigues **✉** !# !# **Date**: 2020-02-17 !# @endnote !# !# @changes !# !# + !# @endchanges !# @bug !# No active bugs reported now !# @endbug !# !# @todo !# ☐
!# @endtodo !# !# @warning !# Now is under CC-GPL License, please see !# © !# @endwarning !# !#--- ---------------------------------------------------------------------------------------- ! character(len=*), intent(in) :: di !# date in long format integer, intent(out) :: yy !# year of date integer, intent(out) :: mm !# month of date integer, intent(out) :: dd !# day of date integer, intent(out) :: hh !# hour of date character(len=4) :: cYear character(len=2) :: cMonth, cDay, cHour integer :: stat !0 1 2 3 !123456789012345678901234567890123 !minutes since 2020-02-13 00:00:00 ! cYear =di(15:18) cMonth=di(20:21) cDay =di(23:24) cHour =di(26:27) !write (*,fmt='(A4,"/",A2,"/",A2," : ",A2)') cYear,cMonth,cDay,cHour read(cYear ,*,iostat=stat) yy read(cMonth,*,iostat=stat) mm read(cDay ,*,iostat=stat) dd read(cHour ,*,iostat=stat) hh !write (*,fmt='(I4.4,"/",I2.2,"/",I2.2," : ",I2.2)') yy,mm,dd,hh hh=hh*100 end subroutine splitDate subroutine adjustMissingValues(icFileType,nVar,nxGrib,nyGrib,nprz_grib2,missingValue & ,levpr_grib2,valuesRead) !# Subroutine to check and fix variable from input !# !# @note !# ![](http://brams.cptec.inpe.br/wp-content/uploads/2015/11/logo-brams.navigation.png "") !# !# **Brief**: This subroutine looks for missing values inside each input variable !# u,v,t,z and r and fix the problem. !# !# If variable is T (Temperature) the subroutine does the replacement using an equation: !# !# \[$T_z=\frac{{P_{z+1}-P_z}}{\rho.C_p}+T_{z+1}$\] !# !# where: !# !# \[$\rho=\frac{P_z}{R.T_{z+1}}$\] !# !# otherwise it just copy the value of variable of level above for the current !# Z position !# !# **Documentation**: !# !# **Author**: Luiz Flavio Rodrigues **✉** !# !# **Date**: 2020-02-19 !# @endnote !# !# @changes !# !# + !# @endchanges !# @bug !# No active bugs reported now !# @endbug !# !# @todo !# ☐
!# @endtodo !# !# @warning !# Now is under CC-GPL License, please see !# © !# @endwarning !# !#--- ---------------------------------------------------------------------------------------- ! include "constants.f90" integer, parameter :: p_tVar=3 integer, intent(in) :: icFileType !# Type of input data integer, intent(in) :: nVar !# Number of incoming variable integer, intent(in) :: nxGrib !# Size of longitudes integer, intent(in) :: nyGrib !# Size of latitudes integer, intent(in) :: nprz_grib2 !# Size of levels real , intent(in) :: levpr_grib2(nprz_grib2) !# Value of level pressure real , intent(in) :: missingValue !# Value for data missing real , intent(inout) :: valuesRead(nxGrib,nyGrib,nprz_grib2) !# Array with variable integer :: i,j,k !# Just counters real :: rho !# $\rho=\frac{P_z}{R.T_{z+1}}$ real :: rhoCpInv !# $rhoCpInv=\frac{1.}{\rho.C_p}$ real :: diffPress !# $diffPress=P_{k}-P_{k+1}$ select case (icFileType) case (3) do k=nprz_grib2-1,1,-1 !Cycle through all levels from highest to surface do j=1,nyGrib !Cycle through all lats and lons do i=1,nxGrib if(valuesRead(i,j,k)==missingValue) then !If var value is missing if(nvar==p_tVar) then !If var is Temperature !Equation rho=levpr_grib2(k)/(c_rgas*valuesRead(i,j,k+1)) rhoCpInv=1./(rho*c_cp) diffPress=levpr_grib2(k)-levpr_grib2(k+1) valuesRead(i,j,k)=valuesRead(i,j,k+1)+diffPress*rhoCpInv else !If var is U,V,Z or R valuesRead(i,j,k)=valuesRead(i,j,k+1) !Just copy the value from level above endif endif enddo enddo enddo case (2) do k=2,nprz_grib2 do j=1,nyGrib do i=1,nxGrib if(valuesRead(i,j,k)==missingValue) then if(nvar==p_tVar) then rho=levpr_grib2(k)/(c_rgas*valuesRead(i,j,k-1)) rhoCpInv=1./(rho*c_cp) diffPress=levpr_grib2(k-1)-levpr_grib2(k) valuesRead(i,j,k)=valuesRead(i,j,k-1)+diffPress*rhoCpInv else valuesRead(i,j,k)=valuesRead(i,j,k-1) endif endif enddo enddo enddo end select end subroutine adjustMissingValues subroutine chem_pressure_stage_grads(n1,n2,n3,nhem,glat,glon,glat2,glon2) use isan_coms use ModDateUtils, only: & date_add_to use chem_isan_coms use chem1_list, only: chemical_mechanism ! intent(in) use mem_chem1, only: CHEM_ASSIM, CHEMISTRY ! intent(in) use rconstants use dump implicit none include "constants.f90" character(len=*),parameter :: sourceName='chem_astp.F90' !Name of source code character(len=*),parameter :: procedureName='**chem_pressure_stage_grads**' !Name of this procedure integer :: n1,n2,n3,nhem,nt real :: glat(n1,n2),glon(n1,n2),glat2(n1,n2),glon2(n1,n2) real :: fnprx,fnpry,grx,gry,gglat,gglon,thmax,thmin & ,xswlon_east,cntlon_east,rr,each_levp_real(maxpr) integer :: i,j,k,lv,ifm,loop,n,iunit,stat real, allocatable :: var(:,:),lat(:,:),lon(:,:) integer :: nx,ny,nz,thisHour integer :: iyyl,imml,iddl,ihhl character(len=256) :: lixo,ctlFileName real :: lonini,dx,latini,dy integer, parameter :: nzmax=256 real :: levs(nzmax) character(len=3) :: cnz integer :: nGradsVars character(len=8),allocatable :: gradsVarsNames(:) character(len=15) :: dataCtl ctlFileName=innpr(1:len_trim(innpr)-3)//'ctl' iErrNumber=dumpMessage(c_tty,c_yes,'','' & ,c_notice,'Reading pressure gridded CTL Grads data '//innpr(1:len_trim(ctlFileName))) open(unit=33,file=trim(ctlFileName) & ,action='read',status='old',form='FORMATTED') read(33,*) lixo read(33,*) lixo read(33,*) lixo read(33,*) lixo,nx,lixo,lonini,dx read(33,*) lixo,ny,lixo,latini,dy read(33,*) lixo,nz,lixo,levs(1:nz) read(33,*) lixo,nt,lixo,dataCtl read(33,*) lixo,nGradsVars allocate(gradsVarsNames(nGradsVars)) do i=1,nGradsVars read(33,*) gradsVarsNames(i) print *,'LFR->',i,gradsVarsNames(i) enddo close(33) allocate(slevs(nz)) do i=1,nz write(slevs(i),fmt='(F5.1)') levs(i) levpr_grib2(i)=levs(i) enddo write(cnz,fmt='(I3.3)') nz write(*,fmt='(A)') '' write(*,fmt='(A)') '¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨' write(*,fmt='("Grads input data inventory for ",A)') DataCtl write(*,fmt='("Nx ......: ",I5.5)') nx write(*,fmt='("Ny ......: ",I5.5)') ny write(*,fmt='("Nz ......: ",I5.5)') nz write(*,fmt='("Init. Lon: ",F8.2,", delta Lon: ",F8.2)') lonini,dx write(*,fmt='("Init. Lat: ",F8.2,", delta Lat: ",F8.2)') latini,dy write(*,fmt='("Final Lon: ",F8.2)') lonini+(nx-1)*dx write(*,fmt='("Final Lat: ",F8.2)') latini+(ny-1)*dy write(*,fmt='("Levels ..: "'//cnz//'(F6.1,1X))') (levpr_grib2(i),i=1,nz) write(*,fmt='(A)') '' write(*,fmt='(A)') '¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨¨' write(*,fmt='("Model limits inventory",A)') '' write(*,fmt='("Nx ......: ",I5.5)') n1 write(*,fmt='("Ny ......: ",I5.5)') n2 write(*,fmt='("Nz ......: ",I5.5)') n3 write(*,fmt='("Init. Lon: ",F8.2,", delta Lon: ",F8.2)') glon(1,1),glon(2,1)-glon(1,1) write(*,fmt='("Init. Lat: ",F8.2,", delta Lat: ",F8.2)') glat(1,1),glat(1,2)-glat(1,1) write(*,fmt='("Final Lon: ",F8.2)') glon(n1,1) write(*,fmt='("Final Lat: ",F8.2)') glat(1,n2) write(*,fmt='(A)') '' if( glon(1,1)lonini+(nx-1)*dx & .or. glat(1,1)latini+(ny-1)*dy) then iErrNumber=dumpMessage(c_tty,c_yes,sourceName,procedureName & ,c_fatal,'Bad domain specification, see limits above!') endif !If the total of levels ir greater than defined in RAMSIN use the defined if(nz>z_max_level) then nprz=z_max_level else nprz=nz endif nprz_grib2=nz nxGrib=nx nyGrib=ny ! ! Increments lat and lon gdatdy=dy gdatdx=dx !Make number of points lon and lat accordingly setup in RAMSIN nprx=nx npry=ny !Setting initial lon and lat xswlon=lonini xswlat=latini inproj=1 if (inproj.eq.1) then ! If necessary, convert longitude specification to values in the range ! [-180.,180.]. xswlon=mod(xswlon+900.,360.)-180. xnelon=mod(xnelon-900.,-360.)+180. if( xswlon<-180.00 .or. xswlon > 180.00 & .or. xnelon<-180.00 .or. xnelon > 180.01 & .or. xswlat< -90.00 .or. xnelat > 90.01) then write(*,fmt='("xswlon: ",F8.2)') xswlon write(*,fmt='("xnelon: ",F8.2)') xnelon write(*,fmt='("xswlat: ",F8.2)') xswlat iErrNumber=dumpMessage(c_tty,c_yes,sourceName,procedureName & ,c_fatal,'Bad domain specification, see info above!') endif fnprx=float(nprx) fnpry=float(npry) ! Set "global domain" flags to determine whether 2 extra rows will be ! added to borders of the input gridded data. iglobew = 1 if doing ! full 360 degrees of longitude. iglobs = 1 if including the south ! pole and if iglobew = 1, and iglobn = 1 if including ! the north pole and if iglobew = 1. iglobew=0 iglobn=0 iglobs=0 if(fnprx*gdatdx > 359.9) then if(abs(xswlon+180.0)<1.e-3) then print*,'When extracting all longitudes of pressure data,' print*,'xswlon must be set to -180.' endif iglobew=1 endif if(iglobew==1 .and. xswlat<-89.9) iglobs=1 if(iglobew==1 .and. xnelat> 89.9) iglobn=1 idatelin=0 if(xnelon.lt.xswlon+.01) idatelin=1 ! Search coarse grid points for any that are outside the bounds of ! the pressure data. If any are found, stop. loop=1 if(nhem>0) loop=2 do ifm = 1,loop do j = 1,n2 do i = 1,n1 if(ifm==1) then gglat=glat(i,j) gglon=glon(i,j) elseif(ifm==2) then gglat=glat2(i,j) gglon=glon2(i,j) endif gry = (gglat - xswlat) / gdatdy + 1. if(gry<2-iglobs .or. gry>fnpry+iglobn-1) then write(*,fmt='("i : ",I4.4)') i write(*,fmt='("j : ",I4.4)') j write(*,fmt='("xswlat : ",F8.2)') xswlat write(*,fmt='("xnelat : ",F8.2)') xnelat write(*,fmt='("gLat : ",F8.2)') gglat write(*,fmt='("iglobs : ", I1)') iglobs write(*,fmt='("iglobn : ", I1)') iglobn write(*,fmt='("fnpry : ",F8.2)') fnpry write(*,fmt='("gdatdy : ",F8.2)') gdatdy write(*,fmt='("gry : ",F8.2)') gry if(gry.lt.2-iglobs) write(*,fmt='(A)') "gry is < 2-iglobs!" if(gry.gt.fnpry+iglobn-1) write(*,fmt='(A)') "gry is > fnpry+iglobn-1!" iErrNumber=dumpMessage(c_tty,c_yes,sourceName,procedureName & ,c_fatal,'Model grid point latitude must be at least 1 ' & //'gridpoint inside pressure data area, see info above!') endif if(idatelin.eq.1.and.gglon.lt.xswlon) gglon=gglon+360. grx=(gglon-xswlon)/gdatdx+1. if(grx < 2-iglobew .or. grx>fnprx+iglobew-1) then write(*,fmt='("i : ",I4.4)') i write(*,fmt='("j : ",I4.4)') j write(*,fmt='("xswlon : ",F8.2)') xswlon write(*,fmt='("xnelon : ",F8.2)') xnelon write(*,fmt='("gLon : ",F8.2)') gglon write(*,fmt='("iglobew: ", I1)') iglobew write(*,fmt='("iglobn : ", I1)') iglobn write(*,fmt='("fnprx : ",F8.2)') fnprx write(*,fmt='("gdatdx : ",F8.2)') gdatdx write(*,fmt='("grx : ",F8.2)') grx if(grx<2-iglobew) write(*,fmt='(A)') "grx is < 2-iglobew!" if(grx>fnprx+iglobew-1) write(*,fmt='(A)') "grx is > fnprx+iglobew-1!" iErrNumber=dumpMessage(c_tty,c_yes,sourceName,procedureName & ,c_fatal,'Model grid point longitude must be at least 1 ' & //'gridpoint inside pressure data area, see info above!') endif enddo enddo enddo endif ! Deallocate memory for the pressure data if(allocated(p_u)) deallocate(p_u) if(allocated(p_v)) deallocate(p_v) if(allocated(p_t)) deallocate(p_t) if(allocated(p_z)) deallocate(p_z) if(allocated(p_r)) deallocate(p_r) if(allocated(p_ur)) deallocate(p_ur) if(allocated(p_vr)) deallocate(p_vr) if(allocated(p_lat)) deallocate(p_lat) if(allocated(p_lon)) deallocate(p_lon) if(allocated(p_slp)) deallocate(p_slp) if(allocated(p_sfp)) deallocate(p_sfp) if(allocated(p_sft)) deallocate(p_sft) if(allocated(p_snow)) deallocate(p_snow) if(allocated(p_sst)) deallocate(p_sst) if(allocated(p_sc)) deallocate(p_sc) if(allocated(p_aer_sc)) deallocate(p_aer_sc) ! Allocate memory for the pressure data allocate(p_u(nprx,npry,nprz)) allocate(p_v(nprx,npry,nprz)) allocate(p_t(nprx,npry,nprz)) allocate(p_z(nprx,npry,nprz)) allocate(p_r(nprx,npry,nprz)) !print *,'LFR: ',nprx,npry,nprz,nspecies_aer_in if(CHEM_ASSIM == 1 .and. nspecies>0 ) allocate(p_sc(nprx,npry,nprz,nspecies)) if(CHEM_ASSIM == 1 .and. nspecies_aer_in>0 ) allocate(p_aer_sc(nprx,npry,nprz,nspecies_aer_in)) ! p_ur,p_vr arrays for rotated winds used later allocate(p_ur(nprx,npry,nprz)) allocate(p_vr(nprx,npry,nprz)) allocate(p_lat(nprx,npry)) allocate(p_lon(nprx,npry)) allocate(p_slp(nprx,npry)) allocate(p_sfp(nprx,npry)) allocate(p_sft(nprx,npry)) allocate(p_snow(nprx,npry)) allocate(p_sst(nprx,npry)) ! Fill with missing in case they are not present p_slp (1:nprx,1:npry)=1.e30 p_sfp (1:nprx,1:npry)=1.e30 p_sft (1:nprx,1:npry)=1.e30 p_snow(1:nprx,1:npry)=1.e30 p_sst (1:nprx,1:npry)=1.e30 iErrNumber=dumpMessage(c_tty,c_yes,'','' & ,c_notice,'Access coarse resolution pressure data') write(*,*) '' ! Call routine to fill pressure arrays from the chosen dataset. call chem_get_press_grads(lonini,dx,latini,dy,nGradsVars,gradsVarsNames) ! Find max-min theta at bottom and top levels thmax=1. thmin=1000. do j=1,npry do i=1,nprx if(p_t(i,j,nprz).lt.1.e20) & thmax=max(thmax,p_t(i,j,nprz)*(p00/pnpr(nprz))**rocp) if(p_t(i,j,1).lt.1.e20) & thmin=min(thmin,p_t(i,j,1)*(p00/pnpr(1))**rocp) enddo enddo print 300,levpr(1),thmin,levpr(nprz),thmax 300 format(//,' Minimum THETA at ',I4,' mb - ',F8.2/ & ,' Maximum THETA at ',I4,' mb - ',F8.2) ! Compute lat-lon at input pressure data points select case (inproj) case (1) do j=1,npry do i=1,nprx p_lat(i,j)=xswlat+(j-1)*gdatdy p_lon(i,j)=xswlon+(i-1)*gdatdx enddo enddo case (2) print*,'lamb-con:',xswlat,xswlon,cntlat,cntlon,gdatdx,gdatdy do j=1,npry do i=1,nprx call lc2ll(cntlat,cntlon,p_lat(i,j),p_lon(i,j) & ,float(i),float(j),xswlat,xswlon,gdatdx) ! print*,i,j,p_lat(i,j),p_lon(i,j) enddo enddo case (3) print*,'polar:',cntlat,cntlon,gdatdx,gdatdy,xswlat,xswlon xswlon_east=xswlon if(xswlon<0.) xswlon_east=xswlon+360. cntlon_east=cntlon if(cntlon<0.) cntlon_east=cntlon+360. do j=1,npry do i=1,nprx call w3fb07(float(i),float(j),xswlat,xswlon_east,gdatdx,cntlon_east & ,p_lat(i,j),p_lon(i,j)) ! print*,i,j,p_lat(i,j),p_lon(i,j) enddo enddo end select end subroutine chem_pressure_stage_grads subroutine chem_get_press_grads(lonini,dx,latini,dy,nGradsVars,gradsVarsNames) use isan_coms use chem_isan_coms #ifdef GRIB2 use wgrib2api #endif use dump, only: & dumpMessage implicit none include "constants.f90" real, intent(in) :: lonini,dx,latini,dy ! logical, external :: checkInside integer, intent(in) :: nGradsVars character(len=8), intent(in) :: gradsVarsNames(nGradsVars) integer :: ithere(maxpr,5),isfthere(5) character*4 field,idat(5) data idat/'T','R','U','V','H'/ character(len=60) :: varName(5) real,allocatable::as(:,:) real, allocatable :: lat(:,:),lon(:,:) real, allocatable :: aux_as(:,:) integer :: i,j,k,nv,nvar,misstot,lv,n,lv_aux integer :: irec,recordLen,iCount,jCount integer :: iAbove(5),iBellow(5) real :: vMax(5),vMin(5) real, allocatable :: dados(:,:,:,:) character(len=256) :: message integer, external :: outRealSize !Allocate array for trimmer area allocate(aux_as(nprx,npry)) iAbove=0 iBellow=0 vMax=uLimit vMin=dlimit allocate(lat(nxGrib,nyGrib),lon(nxGrib,nyGrib),as(nxGrib,nyGrib)) allocate(dados(5+nspecies+nspecies_aer_in,nxGrib,nyGrib,nprz_grib2)) !Making lat and lon arrays do i=1,nxGrib lon(i,:)=lonini+(i-1)*dx !write (88,*) 'Lon:',i,lon(i,1) end do do j=1,nyGrib lat(:,j)=latini+(j-1)*dy !write (89,*) "Lat:",j,lat(1,j) enddo !print *,'INITIAL_LATITUDE: ',INITIAL_LATITUDE,minval(lat) !print *,'FINAL_LATITUDE: ',FINAL_LATITUDE,maxval(lat) !print *,'INITIAL_LONGITUDE:',INITIAL_LONGITUDE,minval(lon) !print *,'FINAL_LONGITUDE: ',FINAL_LONGITUDE,maxval(lon) ! Read upper air fields recordLen=4*nxGrib*nyGrib !srf recordLen=4*nxGrib*nyGrib open(unit=33,file=trim(innpr),& action='read',status='old',form='UNFORMATTED',access='DIRECT', & recl=recordLen) !# reading grads binary and fill variables irec=1 do nvar=1,5+nspecies+nspecies_aer_in do k=1,nprz_grib2 read (33,rec=irec) dados(nvar,1:nxGrib,1:nyGrib,k) irec=irec+1 enddo enddo close(33) do lv=1,nprz_grib2 levpr(lv)=levpr_grib2(lv) lv_aux=lv pnpr(lv_aux)=levpr_grib2(lv)*100 !Do not process levels above z_max_level ! if(lv_aux>z_max_level) cycle do nvar=1,5+nspecies+nspecies_aer_in if(nvar<=5) then aux_as=dados(nvar,:,:,lv)*scale_factor(nvar) else aux_as=dados(nvar,:,:,lv) endif write(message,fmt='("Reading pressure field ",I4," mBar for ",A8)') & int(levpr_grib2(lv)),gradsVarsNames(nvar) iErrNumber=dumpMessage(c_tty,c_yes,'','' & ,c_notice,trim(message)) if(nvar.eq.1) then call prfill(nprx,npry,aux_as,p_u(1,1,lv_aux)) elseif(nvar.eq.2) then call prfill(nprx,npry,aux_as,p_v(1,1,lv_aux)) elseif(nvar.eq.3) then call prfill(nprx,npry,aux_as,p_t(1,1,lv_aux)) elseif(nvar.eq.4) then call prfill(nprx,npry,aux_as,p_z(1,1,lv_aux)) elseif(nvar.eq.5) then call prfill(nprx,npry,aux_as,p_r(1,1,lv_aux)) elseif(nvar > 5 .and. nvar<=5+nspecies) then print *,'LFR-Chem>nvar,ini',nvar,nvar-5 call prfill(nprx,npry,aux_as,p_sc(1,1,lv,nvar-5)) elseif(nvar>5+nspecies) then print *,'LFR-Aer>nvar,ini',nvar,nvar-5-nspecies call prfill(nprx,npry,aux_as,p_aer_sc(1,1,lv,nvar-5-nspecies)) endif enddo enddo if(ccGradsWrite==1) call writeGradsSub4Grads(nGradsVars,gradsVarsNames) !Fill surface fields do nvar=1,5 as=-999.0 select case (nvar) case (1) call prfill(nprx,npry,as,p_slp(1,1)) case (2) call prfill(nprx,npry,as,p_sfp(1,1)) case (3) call prfill(nprx,npry,as,p_sft(1,1)) case (4) call prfill(nprx,npry,as,p_snow(1,1)) case (5) call prfill(nprx,npry,as,p_sst(1,1)) end select enddo ! Check for levels that may be all missing ithere(1:nprz,1:5)=0 isfthere(1:5)=0 do k=1,nprz do j=1,npry do i=1,nprx if(p_t(i,j,k) > 1e20) ithere(k,1)=ithere(k,1)+1 if(p_r(i,j,k) > 1e20) ithere(k,2)=ithere(k,2)+1 if(p_u(i,j,k) > 1e20) ithere(k,3)=ithere(k,3)+1 if(p_v(i,j,k) > 1e20) ithere(k,4)=ithere(k,4)+1 if(p_z(i,j,k) > 1e20) ithere(k,5)=ithere(k,5)+1 enddo enddo enddo do nv=1,5 do k=1,nprz if(ithere(k,nv) < nprx*npry) then ithere(k,nv)=1 else ithere(k,nv)=0 endif enddo enddo misstot=0 do nv=1,5 do k=1,nprz if(ithere(k,nv) == 0) misstot=misstot+1 enddo enddo do j=1,npry do i=1,nprx if(p_slp (i,j) > 1e20) isfthere(1)=isfthere(1)+1 if(p_sfp (i,j) > 1e20) isfthere(2)=isfthere(2)+1 if(p_sft (i,j) > 1e20) isfthere(3)=isfthere(3)+1 if(p_snow(i,j) > 1e20) isfthere(4)=isfthere(4)+1 if(p_sst (i,j) > 1e20) isfthere(5)=isfthere(5)+1 enddo enddo do nv=1,5 if(isfthere(nv) < nprx*npry) then isfthere(nv)=1 else isfthere(nv)=0 endif enddo print*,'------------------------------------------------' print*,' Missing parameters: 0 = all missing' print*,'------------------------------------------------' print '(t20,5(a1,6x))',(idat(n),n=1,5) do k=1,nprz print '(f10.1,t14,5(i7))',pnpr(k),(ithere(k,n),n=1,5) enddo print*,'------------------------------------------------' print*,' Missing surface parameters: 0 = all missing' print*,'------------------------------------------------' print*,' SLP SFP SFT SNOW SST' print '(t14,5(i7))',(isfthere(n),n=1,5) if(misstot.gt.0) then ! Let's see if we can get creative and make up data for the missing fields call press_miss(nprx,npry,nprz,p_u,p_v,p_t,p_z,p_r & ,ithere,maxpr,levpr) ! Check again for missing fields misstot=0 do nv=1,5 do k=1,nprz if(ithere(k,nv).eq.0) misstot=misstot+1 enddo enddo print*,'------------------------------------------------' print*,' After missing parameters check: 0 = all missing' print*,'------------------------------------------------' print '(t20,5(a1,6x))',(idat(n),n=1,5) do k=1,nprz print '(f10.1,t14,5(i7))',pnpr(k),(ithere(k,n),n=1,5) enddo endif deallocate(slevs) !deallocate(mask) !write(*,fmt='(A,I4.4,2(I2.2),I4.4)') 'All done for ',iyear,imonth,idate,ihour return !--(DMK-CCATT-INI)---------------------------------------------------------------- end subroutine chem_get_press_grads subroutine writeGradsSub4Grads(nGradsVars,gradsVarsNames) !Header_vai_aqui! use isan_coms use chem_isan_coms implicit none include "constants.f90" !Parameters (constants) ! Input/Output variables integer, intent(in) :: nGradsVars !,nspecies character(len=8) :: gradsVarsNames(nGradsVars) !Local variables integer :: recordLen, irec, nvar,k character(len=4) :: cYear,cHour character(len=2) :: cDay,cMonth character(len=256) :: gradsFileName character(len=19) :: fileName integer, external :: outRealSize real :: cs(nprx,npry) !Code write(cYear,fmt='(I4.4)') iyear write(cMonth,fmt='(I2.2)') iMonth write(cDay,fmt='(I2.2)') idate write(cHour,fmt='(I4.4)') ihour fileName='ciREad.'//cYear//cMonth//cDay//cHour gradsFileName=trim(ICGRADSPREFIX)//fileName !srf recordLen=4*nprx*npry recordLen=4*nprx*npry open(unit=33,file=trim(gradsFileName)//'.gra',& action='WRITE',status='REPLACE',form='UNFORMATTED',access='DIRECT', & recl=recordLen) !# writing grads binary and fill variables irec=1 do nvar=1,5 do k=1,nprz select case (nvar) case (1) write (33,rec=irec) p_u (1:nprx,1:npry,k) case (2) write (33,rec=irec) p_v (1:nprx,1:npry,k) case (3) write (33,rec=irec) p_t (1:nprx,1:npry,k) case (4) write (33,rec=irec) p_z (1:nprx,1:npry,k) case (5) write (33,rec=irec) p_r (1:nprx,1:npry,k) end select irec=irec+1 enddo enddo do nvar=1,nspecies do k=1,nprz cs=p_sc(:,:,k,nvar) write (33,rec=irec) cs irec=irec+1 enddo enddo do nvar=1,nspecies_aer_in do k=1,nprz cs=p_aer_sc(:,:,k,nvar) write (33,rec=irec) cs irec=irec+1 enddo enddo close(33) open(unit=33,file=trim(gradsFileName)//'.ctl' & ,action='WRITE',status='replace',form='FORMATTED') !writing the name of grads file write(33,*) 'dset ^'//trim(gradsfileName)//'.gra' !writing others infos to ctl write(33,*) 'undef -0.9990000E+34' write(33,*) 'title GRIB NCEP GFS' write(33,*) 'xdef ',nprx,' linear ',xswlon,gdatdx write(33,*) 'ydef ',npry,' linear ',xswlat,gdatdy write(33,*) 'zdef ',nprz,'levels',levpr(1:nprz) write(33,*) 'tdef 1 linear 00:00z01jan2018 1mo' write(33,*) 'vars ',nGradsVars do nvar=1,5+nspecies+nspecies_aer_in write(33,*) gradsVarsNames(nvar),nprz,'99 ',gradsVarsNames(nvar) enddo write(33,*) 'endvars' close(33) end subroutine writeGradsSub4Grads