! This code will read the daily Bootstrap Sea Ice Concentrations from Nimbus-7 ! SMMR and DMSP SSM/I avialable from NSIDC (Comiso, J. 1999) and ! produce intermediate files for WRF-WPS 3.0 every 12 hours for one ! year. !(1) need to change the prefix to point the input file name ! where you have saved the sea-ice concentrations and the number of days ! in Feb if interested in leap year. !(1) Will need to chamnge nx,ny and NREL if working in Northern !Hemisphere !(2) TO COMPILE: pgf90 -Mfree -byteswapio -O seaice-intermd_write.F90 !(3) Move the sea ice files to the directory with the other intermediate !files before running metgrid !(4) You will need to change parameters in the intermediate header for !other seaice data sets integer :: version ! Format version (must =5 for WPS format) ! Use with the north pole grid integer :: nx=275 integer :: ny=205 PARAMETER(JJ=205, II=275, NREL=II*JJ*2) integer*2 data(nx*ny) integer :: i,j,k,DAYS(12),day,dytim integer :: iproj ! Code for projection of data in array: ! 0 = cylindrical equidistant ! 1 = Mercator ! 3 = Lambert conformal conic ! 4 = Gaussian ! 5 = Polar stereographic real :: nlats ! Number of latitudes north of equator ! (for Gaussian grids) real :: xfcst ! Forecast hour of data real :: xlvl ! Vertical level of data in 2-d array real :: startlat, startlon ! Lat/lon of point in array indicated by ! startloc string real :: deltalat, deltalon ! Grid spacing, degrees real :: dx, dy ! Grid spacing, km real :: xlonc ! Standard longitude of projection real :: truelat1, truelat2 ! True latitudes of projection real :: earth_radius ! Earth radius, km real, dimension(275,205) :: slab,slab2 ! The 2-d array holding the data logical :: is_wind_grid_rel ! Flag indicating whether winds are ! relative to source grid (TRUE) or ! relative to earth (FALSE) character (len=8) :: startloc ! Which point in array is given by ! startlat/startlon; set either ! to 'SWCORNER' or 'CENTER ' character (len=9) :: field ! Name of the field character (len=24) :: hdate ! Valid date for data YYYY:MM:DD_HH:00:00 character (len=24),dimension(124) :: fout character (len=25),dimension(31) :: fin character (len=25) :: units ! Units of data character (len=44) :: infile ! Inputfile (binary) with NSIDC sea ice data character (len=2) :: nmon,nmday,hrday !Used with internal write for file names character (len=32) :: map_source ! Source model / originating center character (len=46) :: desc ! Short description of data character (len=14) :: sfix !file name suffix character (len=26) :: prefix !input file name prefix character (len=26) :: opfix !output file name prefix character (len=35) :: outfile character (len=1) :: sp,ul DATA DAYS/31,28,31,30,31,30,31,31,30,31,30,31/ DATA& startloc,startlat,startlon,dx,dy,xlonc,truelat1,earth_radius& /'CENTER ', 80.50, 160.00,50,50, -114.0, 80.5, 6370.000/ DATA field,units,map_source,desc,iproj,version,xlvl,xfcst & /'SEAICE','Fraction (%)','NSIDC Comiso, J. 1999',& 'Bootstrap Sea Ice Concentrations Nimbus-7 SMMR and DMSP SSM/I',& 5, 5, 200100,0/ ! Prefix used to build input file name based on location of files ! downloaded fron NSIDC; opfix is the output prefix; cuurently set ! to produce intermediate files for 1998 prefix='../DATA/NSIDC/2007/bt_2007' sfix='_f13_v02_n.bin' opfix='../DATA/INTMDT/SEAICE:2007' sp ='-' ul ='_' print*,'startloc, startlat, startlon',startloc,startlat,startlon print*,'xlonc,truelat1,truelat2, dx, dy ',xlonc,truelat1,dx,dy print*,'Proj, Version, earth_radius',iproj,version,earth_radius DO mon=1,12 !Changes this to DO mon=1,12 for full year if(mon.LT.10)then write(nmon,3) mon else write(nmon,4) mon end if do day=1,DAYS(mon) if(day.LT.10)then write(nmday,3) day else write(nmday,4) day end if 3 format('0',I1) 4 format(I2) infile=prefix//nmon//nmday//sfix print*," INFILE ", infile OPEN(50,FILE=infile,status='old',CONVERT='LITTLE_ENDIAN'& ,recl=nrel, access='DIRECT',form='unformatted') read(50,rec=1) data close(50) do 100 i=1,ii do 100 j=1,jj L=ii*(j-1)+i IF (data(L).LT.1050) THEN x=float(data(L)) slab(i,j)=x/10. ! print *, x, slab else IF (data(L).EQ.1100) THEN slab(i,j)=100 else slab(i,j)=-999. ENDIF ENDIF 100 continue write(66,77) ((slab(i,j),j=1,jj),i=1,ii) 77 format(80(f8.2,',')) dytim=0 DO n=1,2 if(dytim.LT.10)then write(hrday,3) dytim else write(hrday,4) dytim end if dytim=dytim+12 outfile=opfix//sp//nmon//sp//nmday//ul//hrday open(unit=11,file=outfile,form='unformatted') ! for output SEAICE files print7,"IN ",infile," OUT ",outfile write(hdate,39) outfile(6:20),":00:00" 39 format(A13,A6) 7 format(A3,A55,A5,A24) ! 1) WRITE FORMAT VERSION write(unit=11) version ! 2) WRITE METADATA ! Polar stereographic write(unit=11) hdate, xfcst, map_source, field, & & units, desc, xlvl, nx, ny, iproj write(unit=11) startloc, startlat, startlon, dx, dy, & & xlonc, truelat1, earth_radius print*,hdate, xfcst, map_source, field,units print*,startloc, startlat, startlon, dx, dy,nx,ny ! 3) WRITE WIND ROTATION FLAG write(unit=11) is_wind_grid_rel ! 4) WRITE 2-D ARRAY OF DATA do i=1,II k=1 do j=JJ,1,-1 slab2(i,k)=slab(i,j) k=k+1 end do end do write(unit=11) slab2 close(11) END DO ! End of day end do ! End of days in the month END DO ! End of months in the year end