program fmain !------------------------------------------------------------------------------- ! PURPOSE: ! o given a SCRIP map matrix data file, map a field ! ! NOTES: ! o all output data is base on the "_a" grid, the "_b" grid is ignored !------------------------------------------------------------------------------- implicit none include 'netcdf.inc' integer,parameter :: R8 = selected_real_kind(12) ! 8 byte real integer,parameter :: R4 = selected_real_kind( 6) ! 4 byte real integer,parameter :: RN = kind(1.0) ! native real integer,parameter :: I8 = selected_int_kind (13) ! 8 byte integer integer,parameter :: I4 = selected_int_kind ( 6) ! 4 byte integer integer,parameter :: IN = kind(1) ! native integer integer,parameter :: CS = 80 ! short char integer,parameter :: CL = 256 ! long char integer,parameter :: CX = 512 ! extra-long char integer :: n ! index integer :: nargs ! number of arguments integer, external :: iargc ! number of arguments function character(LEN=512) :: arg ! input argument character(LEN=512) :: cmdline ! input command line character(LEN=512) :: fmap ! file name ( input nc file) character(LEN=512) :: fn_out,fn_in ! temporary character(LEN=512) :: var_out,var_in ! temporary character(LEN=512) :: usercomment ! user comment character(LEN= 8) :: cdate ! wall clock date character(LEN=10) :: ctime ! wall clock time !---------------------------------------------------- fmap = 'null' fn_in = 'null' var_in = 'null' fn_out = 'null' var_out = 'null' usercomment = 'null' nargs = iargc() if (nargs == 0) then write(6,*)'invoke map_field -h for usage' stop end if cmdline = 'map_field ' n = 1 do while (n <= nargs) arg = ' ' call getarg (n, arg) n = n + 1 select case (arg) case ('-m') ! input mapping file call getarg (n, arg) n = n + 1 fmap = trim(arg) cmdline = trim(cmdline) // ' -m ' // trim(arg) case ('-of') ! output filename call getarg (n, arg) n = n + 1 fn_out = trim(arg) cmdline = trim(cmdline) // ' -of ' // trim(arg) case ('-ov') ! output variable name call getarg (n, arg) n = n + 1 var_out = trim(arg) cmdline = trim(cmdline) // ' -ov ' // trim(arg) case ('-if') ! input filename call getarg (n, arg) n = n + 1 fn_in = trim(arg) cmdline = trim(cmdline) // ' -if ' // trim(arg) case ('-iv') ! input variable name call getarg (n, arg) n = n + 1 var_in = trim(arg) cmdline = trim(cmdline) // ' -iv ' // trim(arg) case ('-c') ! user comment call getarg (n, arg) n = n + 1 usercomment = trim(arg) case ('-h') call usage_exit (' ') case default write (6,*) 'Argument ', arg,' is not known' call usage_exit (' ') end select end do if (fmap == 'null' .or. fn_out == 'null' .or. fn_in== 'null') then call usage_exit ('Must specify all the following arguments') end if call date_and_time(cdate,ctime) call map_field (fmap, fn_in, var_in, fn_out, var_out, usercomment) !-------------- contains !-------------- subroutine map_field(fmap, fn_in, var_in, fn_out, var_out, usercomment) implicit none !--- includes --- include 'netcdf.inc' ! netCDF defs character(LEN=*), intent(in) :: fmap ! file name ( input nc file) character(LEN=*), intent(in) :: fn_in ! file name character(LEN=*), intent(in) :: var_in ! var name character(LEN=*), intent(in) :: fn_out ! file name character(LEN=*), intent(in) :: var_out ! var name character(LEN=*), intent(in) :: usercomment ! user comment from namelist !--- domain data --- integer :: nb ! size of 1d domain integer :: nbi ! size of i-axis of 2d domain integer :: nbj ! size of j-axis of 2d domain integer :: na ! size of source array integer :: nai ! size of i-axis of 2d domain integer :: naj ! size of j-axis of 2d domain integer :: nd ! size of fld1 real(r8) ,pointer :: xc( :) ! x-coords of center real(r8) ,pointer :: yc( :) ! y-coords of center real(r8) ,pointer :: xv(:,:) ! x-coords of verticies real(r8) ,pointer :: yv(:,:) ! y-coords of verticies real(r8) ,pointer :: area(:) ! cell area real(r8) ,pointer :: fld1(:) ! fld1 real(r8) ,pointer :: fld12(:,:) ! fld1 2d real(r8) ,pointer :: fld2(:) ! fld2 integer ,pointer :: src_grid_dims(:) integer ,pointer :: dst_grid_dims(:) !--- for mapping --- integer :: ns ! size of wgts list integer ,pointer :: col ( :) ! column index integer ,pointer :: row ( :) ! row index real(r8),pointer :: S ( :) ! wgts !--- local --- character(LEN=CL) :: flong ! long name character(LEN=CL) :: str_grido ! global attribute str - grid_file_ocn character(LEN=CL) :: str_grida ! global attribute str - grid_file_atm integer :: fid ! nc file ID integer :: i,j,k ! generic indicies integer :: attnum ! attribute number character(LEN=CL) :: units ! netCDF attribute name string integer :: vid ! nc variable ID integer :: did ! nc dimension ID integer :: dimids(2) ! nc dimension ID integer :: rcode ! routine return error code integer :: dst_grid_rank, src_grid_rank real(r8),parameter :: pi = 3.14159265358979323846 real(r8),parameter :: c0 = 0.00000000000000000000 real(r8),parameter :: c1 = 1.00000000000000000000 real(r8),parameter :: c90 = 90.0000000000000000000 !--- formats --- ! character(LEN=*),parameter :: F00 = "(120a )" ! character(LEN=*),parameter :: F02 = "(a,5i6,i12)" ! character(LEN=*),parameter :: F10=& ! & "('Data created: 'i4,'-',i2,2('-',i2),' ',i2,2('-',i2),' ')" !---------------------------------------------------------------------------- write(6,*) 'fmap = ',trim(fmap) write(6,*) 'fn_in = ',trim(fn_in) write(6,*) 'var_in = ',trim(var_in) write(6,*) 'fn_out = ',trim(fn_out) write(6,*) 'var_out= ',trim(var_out) write(6,*) 'usercomment= ',trim(usercomment) !---------------------------------------------------------------------------- write(6,*) ' ' write(6,*) 'input SCRIP data...' !---------------------------------------------------------------------------- write(6,*) ' ' write(6,*) 'input file = ',fmap(1:len_trim(fmap)) call check_ret(nf_open(fmap(1:len_trim(fmap)),NF_NOWRITE,fid)) write(6,*) 'open ',trim(fmap) str_grido = 'unknown' str_grida = 'unknown' rcode = nf_get_att_text(fid, NF_GLOBAL, 'grid_file_ocn', str_grido) if ( rcode == nf_enotatt ) then rcode = nf_get_att_text(fid, NF_GLOBAL, 'grid_file_src', str_grido) end if if ( rcode == nf_enotatt ) then rcode = nf_get_att_text(fid, NF_GLOBAL, 'domain_a', str_grido) end if rcode = nf_get_att_text(fid, NF_GLOBAL, 'grid_file_atm', str_grida) if ( rcode == nf_enotatt ) then rcode = nf_get_att_text(fid, NF_GLOBAL, 'grid_file_dst', str_grida) end if if ( rcode == nf_enotatt ) then rcode = nf_get_att_text(fid, NF_GLOBAL, 'domain_b', str_grida) end if write(6,*) 'grid_file_ocn= ',trim(str_grido) write(6,*) 'grid_file_atm= ',trim(str_grida) !---------------------------------------------- ! get domain info !---------------------------------------------- call check_ret(nf_inq_dimid (fid, 'n_b', did)) call check_ret(nf_inq_dimlen(fid, did , nb)) rcode = nf_inq_dimid (fid, 'ni_b', did) if (rcode == 0) then call check_ret(nf_inq_dimlen(fid, did, nbi)) else nbi = n end if rcode = nf_inq_dimid (fid, 'nj_b', did) if (rcode == 0) then call check_ret(nf_inq_dimlen(fid, did, nbj)) else nbj = 1 end if if (nbi == 1 .and. nbj == 0) then nbi = n nbj = 1 end if call check_ret(nf_inq_dimid (fid, 'n_s', did)) call check_ret(nf_inq_dimlen(fid, did , ns)) call check_ret(nf_inq_dimid (fid, 'n_a', did)) call check_ret(nf_inq_dimlen(fid, did , na)) call check_ret(nf_inq_dimid (fid, 'dst_grid_rank', did)) call check_ret(nf_inq_dimlen(fid, did, dst_grid_rank)) call check_ret(nf_inq_dimid (fid, 'src_grid_rank', did)) call check_ret(nf_inq_dimlen(fid, did, src_grid_rank)) allocate(src_grid_dims(src_grid_rank), dst_grid_dims(dst_grid_rank)) call check_ret(nf_inq_varid(fid,'src_grid_dims', vid)) call check_ret(nf_get_var_int(fid,vid,src_grid_dims)) call check_ret(nf_inq_varid(fid,'dst_grid_dims', vid)) call check_ret(nf_get_var_int(fid,vid,dst_grid_dims)) if (dst_grid_rank == 2) then nbi = dst_grid_dims(1) nbj = dst_grid_dims(2) else nbi = dst_grid_dims(1) nbj = 1 end if nb = nbi * nbj if (src_grid_rank == 2) then nai = src_grid_dims(1) naj = src_grid_dims(2) else nai = src_grid_dims(1) naj = 1 end if na = nai * naj write(6,*)'na,nai,naj,nb,nbi,nbj,ns= ',na,nai,naj,nb,nbi,nbj,ns allocate(fld1(na)) allocate(fld2(nb)) allocate(col(ns)) allocate(row(ns)) allocate(S(ns)) !--- read weights --- call check_ret(nf_inq_varid(fid,'col', vid )) call check_ret(nf_get_var_int(fid,vid, col )) call check_ret(nf_inq_varid(fid,'row', vid )) call check_ret(nf_get_var_int(fid,vid, row )) call check_ret(nf_inq_varid(fid,'S', vid )) call check_ret(nf_get_var_double(fid,vid,S)) call check_ret(nf_close(fid)) !--- read fld1 --- write(6,*) ' ' write(6,*) 'input file = ',fn_in(1:len_trim(fn_in)) call check_ret(nf_open(fn_in(1:len_trim(fn_in)),NF_NOWRITE,fid)) write(6,*) 'open ',trim(fn_in) call check_ret(nf_inq_varid(fid,trim(var_in), vid )) call check_ret(nf_inq_varndims(fid,vid,nd)) if (nd == 1) then call check_ret(nf_inq_vardimid(fid,vid,dimids)) call check_ret(nf_inq_dimlen(fid,dimids(1),nai)) if (nai /= na) then write(6,*) 'error nai size ',nai,na stop endif call check_ret(nf_get_var_double(fid,vid,fld1 )) elseif (nd == 2) then call check_ret(nf_inq_vardimid(fid,vid,dimids)) call check_ret(nf_inq_dimlen(fid,dimids(1),nai)) call check_ret(nf_inq_dimlen(fid,dimids(2),naj)) if (nai*naj /= na) then write(6,*) 'error nai naj size ',nai,naj,nai*naj,na stop endif allocate(fld12(nai,naj)) call check_ret(nf_get_var_double(fid,vid,fld12 )) n = 0 do j = 1,naj do i = 1,nai n = n + 1 fld1(n) = fld12(i,j) enddo enddo deallocate(fld12) else write(6,*) 'error nd ',nd stop endif call check_ret(nf_close(fid)) !---------------------------------------------------------------------------- write(6,*) 'compute fld2' !---------------------------------------------------------------------------- fld2 = c0 do k = 1,ns fld2(row(k)) = fld2(row(k)) + fld1(col(k))*S(k) enddo !----------------------------------------------------------------- ! create a new nc files !----------------------------------------------------------------- write(6,*) ' ' write(6,*) 'output output file...' units = 'unitless' flong = 'mapped '//trim(var_in) call check_ret(nf_create(fn_out(1:len_trim(fn_out)),NF_CLOBBER,fid)) write(6,*) 'write ',trim(fn_out) call write_file(fid, fn_out, units, nb, nbi, nbj, & fld2, flong, fmap, str_grido, str_grida) call check_ret(nf_close(fid)) write(6,*) 'successfully created domain file ', trim(fn_out) end subroutine map_field !=========================================================================== subroutine check_ret(ret) ! Check return status from netcdf call implicit none integer, intent(in) :: ret if (ret /= NF_NOERR) then write(6,*)'netcdf error with rcode = ', ret,' error = ', nf_strerror(ret) call abort() end if end subroutine check_ret subroutine usage_exit (arg) implicit none character(len=*) :: arg if (arg /= ' ') then write (6,*) arg end if write(6,*) ' Purpose:' write(6,*) ' Map data from one grid to another given a mapping file ' write(6,*) ' and a source field in another netcdf file ' write(6,*) ' ' write(6,*) ' Usage: ' write(6,*) ' map_field -m ' write(6,*) ' -if ' write(6,*) ' -iv ' write(6,*) ' -of ' write(6,*) ' -ov ' write(6,*) ' [-c ]' write(6,*) ' ' write(6,*) ' Where: ' write(6,*) ' filemap = input mapping filename' write(6,*) ' input_file = input file where input field is read' write(6,*) ' input_varname = name of the variable in the input_file to read' write(6,*) ' output_file = where mapped field is written' write(6,*) ' input_varname = name of the mapped field variable in the output_file' write(6,*) ' usercomment = optional, netcdf global attribute (character string)' write(6,*) ' ' write(6,*) ' NOTE that the output_file is always clobbered when using this tool' write(6,*) ' ' stop end subroutine usage_exit !=========================================================================== subroutine write_file(fid, fout, units, n, ni, nj, & fld2, fld2long, fmap, str_grido, str_grida) implicit none !--- includes --- include 'netcdf.inc' ! netCDF defs integer , intent(in) :: fid ! nc file ID character(LEN=*), intent(in) :: fout ! file name ( output nc file) character(LEN=*), intent(in) :: units ! netCDF attribute name string integer , intent(in) :: n ! size of 1d domain integer , intent(in) :: ni ! size of i-axis of 2d domain integer , intent(in) :: nj ! size of j-axis of 2d domain real(r8) , pointer :: fld2(:) ! output field character(LEN=*), intent(in) :: fld2long ! global attribute str - grid_file_ocn character(LEN=*), intent(in) :: fmap ! global attribute str - grid_file_ocn character(LEN=*), intent(in) :: str_grido ! global attribute str - grid_file_ocn character(LEN=*), intent(in) :: str_grida ! global attribute str - grid_file_atm !--- local --- character(LEN=CL) :: host ! hostname of machine running on character(LEN=CL) :: str ! fixed length char string character(LEN=CL) :: str_title ! global attribute str - title character(LEN=CL) :: str_source ! global attribute str - source character(LEN=CL) :: user ! user name integer :: strlen ! (trimmed) length of string integer :: vid ! nc variable ID integer :: did ! nc dimension ID integer :: vdid(3) ! vector of nc dimension ID integer :: rcode ! routine return error code character(*),parameter:: version = 'SVN $Id: map_field.F90 39483 2012-08-16 17:35:52Z mlevy@ucar.edu $' ! global attributes str = 'map_field data: ' call check_ret(nf_put_att_text(fid,NF_GLOBAL,'title' ,len_trim(str),str)) str = 'CF-1.0' call check_ret(nf_put_att_text(fid,NF_GLOBAL,'Conventions',len_trim(str),str)) str = trim(version) call check_ret(nf_put_att_text(fid,NF_GLOBAL,'source_code',len_trim(str),str)) str = ' $URL: https://svn-ccsm-models.cgd.ucar.edu/tools/mapping/map_field/trunk/src/map_field.F90 $' call check_ret(nf_put_att_text(fid,NF_GLOBAL,'SVN_url',len_trim(str),str)) #ifdef OPT str = 'TRUE' #else str = 'FALSE' #endif call check_ret(nf_put_att_text(fid,NF_GLOBAL,'Compiler_Optimized',len_trim(str),str)) call sys_getenv('HOST',host,rcode) if (rcode == 0) then call check_ret(nf_put_att_text(fid,NF_GLOBAL,'hostname' ,len_trim(host),host)) else call sys_getenv('HOSTNAME',host,rcode) if (rcode == 0) then call check_ret(nf_put_att_text(fid,NF_GLOBAL,'hostname' ,len_trim(host),host)) else write(6,*) 'WARNING: could not determine hostname, so that information will not be stored in netCDF attribute. To avoid this warning in the future, set environment variable HOST or HOSTNAME.' end if end if call sys_getenv('LOGNAME',user,rcode) if (rcode /= 0) then write(6,*) ' ERROR: getting LOGNAME' stop end if str = 'created by '//trim(user)//', '//cdate(1:4)//'-'//cdate(5:6)//'-'//cdate(7:8) & & //' '//ctime(1:2)//':'//ctime(3:4)//':'//ctime(5:6) call check_ret(nf_put_att_text(fid,NF_GLOBAL,'history' ,len_trim(str),str)) str = trim(fout) call check_ret(nf_put_att_text(fid,NF_GLOBAL,'source' ,len_trim(str),str)) str = trim(fmap) call check_ret(nf_put_att_text(fid,NF_GLOBAL,'map_file',len_trim(str),str)) str = trim(str_grido) call check_ret(nf_put_att_text(fid,NF_GLOBAL,'map_grid_file_ocn',len_trim(str),str)) str = trim(str_grida) call check_ret(nf_put_att_text(fid,NF_GLOBAL,'map_grid_file_atm',len_trim(str),str)) str = trim(usercomment) if ( str(1:4) /= 'null' ) then call check_ret(nf_put_att_text(fid,NF_GLOBAL,'user_comment',len_trim(str),str)) end if !----------------------------------------------------------------- ! dimension data !----------------------------------------------------------------- call check_ret(nf_def_dim(fid, 'n' , n , did)) ! # of points total call check_ret(nf_def_dim(fid, 'ni', ni, did)) ! # of points wrt i vdid(1) = did call check_ret(nf_def_dim(fid, 'nj', nj, did)) ! # of points wrt j vdid(2) = did !----------------------------------------------------------------- ! define data -- coordinates, input grid !----------------------------------------------------------------- call check_ret(nf_def_var (fid,trim(var_out),NF_DOUBLE ,2,vdid,vid)) str = trim(fld2long) call check_ret(nf_put_att_text(fid,vid,"long_name",len_trim(str),str)) str = trim(units) call check_ret(nf_put_att_text(fid,vid,"units" ,len_trim(str),str)) call check_ret(nf_enddef(fid)) call check_ret(nf_inq_varid(fid,trim(var_out),vid)) call check_ret(nf_put_var_double(fid, vid ,fld2)) end subroutine write_file SUBROUTINE sys_getenv(name, val, rcode) IMPLICIT none !----- arguments ----- character(*) ,intent(in) :: name ! env var name character(*) ,intent(out) :: val ! env var value integer,intent(out) :: rcode ! return code !----- local ----- integer :: lenname ! length of env var name integer :: lenval ! length of env var value character(len=80) :: tmpval ! temporary env var value !----- formats ----- character(*),parameter :: subName = '(shr_sys_getenv) ' character(*),parameter :: F00 = "('(shr_sys_getenv) ',4a)" !------------------------------------------------------------------------------- ! PURPOSE: an architecture independant system call !------------------------------------------------------------------------------- lenname=len_trim(name) #if (defined IRIX64 || defined CRAY || defined UNICOSMP) call pxfgetenv(name, lenname, val, lenval, rcode) #elif (defined AIX || defined OSF1 || defined SUNOS || defined LINUX || defined NEC_SX) #ifdef F2003 call GET_ENVIRONMENT_VARIABLE(trim(name),value=val, length=lenval, status=rcode) #else call getenv(trim(name),tmpval) val=trim(tmpval) rcode = 0 if (len_trim(val) == 0 ) rcode = 1 if (len_trim(val) > CL) rcode = 2 #endif #else write(*,F00) 'ERROR: no implementation of getenv for this architecture' stop subname//'no implementation of getenv for this machine' #endif END SUBROUTINE sys_getenv end program fmain