IDL CMIP6 NC格式数据处理

 下载了一些模型包括以下：

AWI-CM-1-1-MR

ACCESS-CM2

BCC-CSM2-MR

CESM2

CESM2-WACCM

CMCC-CM2-SR5

CMCC-ESM2

CNRM-CM6-1

CNRM-ESM2-1

CanESM5

EC-Earth3

EC-Earth3-Veg

EC-Earth3-Veg-LR

IPSL-CM6A-LR

MIROC-ES2L

MIROC6

MPI-ESM1-2-HR

MPI-ESM1-2-LR

MRI-ESM2-0

NorESM2-LM

NorESM2-MM

EC-Earth3-CC

KIOST-ESM

我自己下载的是siconc和sithick两个变量，但是本程序已经可以根据文件名自动判别变量名称了，而且存在SIMon，和SIday两种frequency的。

所以要求程序能够自动处理不同模型、变量不同、频率不同的情况。

程序只把数据提取成envi标准格式，并且也存了envi的hdr，文件名需要保留除日期之外的全部，再加上日期/月份，存储的文件是月/日的。

最后也保存了lon和lat

我觉得比较难做的主要是，各种日期的识别，因为有考虑闰年的，不考虑闰年的之类，因此都单独用了函数来计算。

pro read_CMIP6_model_data
  infilepath = 'H:\mission\LNG\Model\MPI-ESM1-2-LR\ssp585\daily\sicon'
  outfilepath = 'H:\mission\LNG\Model\MPI-ESM1-2-LR\ssp585\daily\sicon\hdr\'
  nv = envi(/HEADLESS)

  CD,infilepath
  thesefiles = FILE_SEARCH('*.nc',count = numFiles)
  print,'number of files found:',numFiles
  ;处理数据

  FOR fidx=0, numFiles-1 DO BEGIN; ;不要忘记修改这个
    filename=file_basename(thesefiles[fidx],'.nc')
    str=STRSPLIT(filename, '_', /EXTRACT) ;siconc_SIday_AWI-CM-1-1-MR_historical_r1i1p1f1_gn_18500101-18501231
    ;获取的是一期的数据
    file_ID = NCDF_OPEN(theseFiles[fidx],/NOWRITE)   ;open netCDF file for READ only
    Tag = NCDF_INQUIRE(file_ID)
    timeid = NCDF_VARID(file_ID,'time')             ;initial_time0initial_time0_hours  initial_time0
    NCDF_VARGET, file_ID, timeid, time              ;获取 时间维变量
    nt = N_ELEMENTS(time)
    NCDF_ATTGET,file_ID, timeid,'calendar',timeinfo
    timeinfo=string(timeinfo)
    NCDF_ATTGET,file_ID, timeid,'units',timeunits
    timeunits=string(timeunits)
    timeunits=STRSPLIT(timeunits, ' ', /EXTRACT)
    
    ;print,time
    latid = NCDF_VARID(file_ID,'lat')               ;g0_lat_2
    if latid eq -1 then begin
      latid = NCDF_VARID(file_ID,'latitude')  
      if latid eq -1 then begin
        latid = NCDF_VARID(file_ID,'nav_lat')  
      endif
    endif
    NCDF_VARGET, file_ID, latid, latitude           ;获取 维度变量
    
    nlat = N_ELEMENTS(latitude)                     ;纬向维长度
    ;print,'dsklfjdsjfkldjsfkljl'+latitude
    lonid = NCDF_VARID(file_ID,'lon') 
    if lonid eq -1 then begin
      lonid = NCDF_VARID(file_ID,'longitude')
      if lonid eq -1 then begin
        lonid = NCDF_VARID(file_ID,'nav_lon')
      endif
    endif               ;g0_lon_3
    NCDF_VARGET, file_ID, lonid, longitude           ;获取 经度变量
    nlon = N_ELEMENTS(longitude)                     ;径向维长度

    if size(latitude,/N_DIMENSIONS) eq 1 then begin ;如果经纬度为一维的，就要建立成二维的
      templa=fltarr(nlon,nlat)
      templo=fltarr(nlon,nlat)
;      templa=MAKE_ARRAY(nlat,nlon,/FLOAT,VALUE =latitude[*])
      for i=0,nlat-1 do begin
        templa[*,i]=latitude[i]
      endfor
      for i=0,nlon-1 do begin
        templo[i,*]=longitude[i]
      endfor
      
      latitude=templa
      longitude=templo
    endif
    
    var_ID = NCDF_VARID(file_ID,str[0] )  ;获取变量IDvar_result.Name  siconc'sithick'
    NCDF_VARGET, file_ID, Var_ID, var_Data        ;获取变量数据
    
    NCDF_ATTGET,file_ID, var_ID,'missing_value',varinfo
    varinfo=float(varinfo)

    timeunits=STRSPLIT(timeunits[2], '-', /EXTRACT)
    year=long(timeunits[0]);1850;long(STRMID(str[6],0,4))
    month=long(timeunits[1]);long(STRMID(str[6],4,2))
    day=long(timeunits[2])

    
    NCDF_CLOSE, file_ID
    print,'  All parameters of file ' + thesefiles[fidx] + ' have been writen out!'
    index=where(latitude[*,*] ge 47,count) ;我只要纬度>47的数据
    ind = ARRAY_INDICES(latitude,index)

    ;    输出经纬度
    dim=size(latitude[index],/DIMENSIONS)
    if N_elements(dim) eq 1 then begin
      dim=[dim,1]
    endif
    envi_setup_head,FNAME=infilepath+'\lat',NS=dim[0] ,NL=dim[1],NB=1,DATA_TYPE=size(latitude,/type),INTERLEAVE=0,/OPEN,/WRITE ;5为double
    openw,lunou,infilepath+'\lat',/get_lun
    writeu,lunou,latitude[index]
    close,lunou
    free_lun,lunou
    envi_setup_head,FNAME=infilepath+'\lon',NS=dim[0] ,NL=dim[1],NB=1,DATA_TYPE=size(longitude,/type),INTERLEAVE=0,/OPEN,/WRITE
    openw,lunou,infilepath+'\lon',/get_lun
    writeu,lunou,longitude[index]
    close,lunou
    free_lun,lunou


;通过时间变量的计算判定日期
    for x=0,nt-1 do begin
      if timeinfo eq '360_day' then begin
        tempresult=calendar_360(year,month,day,time[x])
        yr=tempresult[0]
        mon=tempresult[1]
        dy=tempresult[2]
      endif else if timeinfo eq 'noleap' then begin
        tempresult=calendar_noleap(year,month,day,time[x])
        yr=tempresult[0]
        mon=tempresult[1]
        dy=tempresult[2]
      endif else begin
        caldat,julday(month,day,year)+time[x]-0.5,mon,dy,yr;这个是proleptic_gregorian或者儒略历(Julian) thick需要在filename后面添加这个
;        print,mon,dy,yr
;        gregorian and julday For dates between 1 Jan CE and 15 Oct 1582 the two calendar systems differ by up to 10 days. For dates on or after 15 Oct 1582 the two calendar systems are identical.
      endelse
      
      
      if yr lt 1980 then begin
        continue
      endif
      
      outfilename = outfilepath +String(yr,FORMAT='(I04)')+'\'+strjoin(str[0:5],'_')+'_'+String(yr,FORMAT='(I04)')+String(mon,FORMAT='(I02)');+String(dy,FORMAT='(I02)');STRCOMPRESS(outfilepath +  'smos_sea_ice_thickness_'+year+'_'+month+'_'+day+'.tif',/remove_all)
      ;存储数据
      if str[1] eq 'SIday' then begin
        outfilename = outfilename+String(dy,FORMAT='(I02)')
        if mon eq 2 and dy eq 29 then begin ;不要闰年的2月29号数据
          continue
        endif
        if yr gt 2059 then begin
          continue
        endif
      endif
      
;      temp=[]
      if size(var_Data,/N_DIMENSIONS) eq 3 then begin ;这是由于变量不是像AWI，而是两维+1维时间
        dim=size(var_Data,/DIMENSIONS)
        temp=fltarr(dim[0],dim[1])
        temp[*,*]=var_Data[*,*,x]
        temp=temp[index];var_Data[ind[1],ind[0],x]
      endif else begin
        temp=var_Data[index,x]
      endelse
      miss=where(temp eq varinfo,count)
      if count gt 0 then begin
        temp[miss]=-1
      endif
;      temp1=var_Data[ind[0],ind[1],x]
;      print,max(temp),mean(temp),var_Data[1234,137,x];,temp[100,37]
      dim=size(temp,/DIMENSIONS)
      if N_elements(dim) eq 1 then begin
        dim=[dim,1]
      endif
      envi_setup_head,FNAME=outfilename,NS=dim[0] ,NL=dim[1],NB=1,DATA_TYPE=size(temp,/type),INTERLEAVE=0,/OPEN,/WRITE ;4为float
      openw,lunou,outfilename,/get_lun
      writeu,lunou,temp;
      close,lunou
      free_lun,lunou
      ;消除内存
      temp=!Null
      delvar,temp
      print,outfilename
      

      
    endfor
    ;    输出经纬度
    dim=size(latitude[index],/DIMENSIONS)
    if N_elements(dim) eq 1 then begin
      dim=[dim,1]
    endif
    envi_setup_head,FNAME=infilepath+'\lat',NS=dim[0] ,NL=dim[1],NB=1,DATA_TYPE=size(latitude,/type),INTERLEAVE=0,/OPEN,/WRITE ;5为double
    openw,lunou,infilepath+'\lat',/get_lun
    writeu,lunou,latitude[index]
    close,lunou
    free_lun,lunou
    envi_setup_head,FNAME=infilepath+'\lon',NS=dim[0] ,NL=dim[1],NB=1,DATA_TYPE=size(longitude,/type),INTERLEAVE=0,/OPEN,/WRITE
    openw,lunou,infilepath+'\lon',/get_lun
    writeu,lunou,longitude[index]
    close,lunou
    free_lun,lunou
    ;    var_Data[where(~finite(var_Data))] = -1
    var_Data=!Null
    delvar,var_Data
    
    

    print,thesefiles[fidx],' is finished'
  ENDFOR


  nv.close


end
function calendar_360,yr,mon,dy,time
  tdy=dy+time mod 30
  tmon=long(time/30)
  tyr=yr+long(tmon/12)
  tmon=mon+tmon mod 12
  return,[tyr,tmon,tdy]
end
FUNCTION calendar_noleap,yr,mon,dy,time
  ;不考虑闰年
  ;
  months_days = [31,28,31,30,31,30,31,31,30,31,30,31]
  time=time+0.5
  ;  IF ISLEAPYEAR(year) EQ 1 THEN months_days[1]=28;29 ;注意我这边由于闰年的28都被排除了，所以DOY不计算29的情况
  tyr=yr+long(time/365)
  tdy=time mod 365
  for month=0,12-1 do begin
    if tdy le total(months_days[0:month]) then begin
      break
    endif
  endfor
  tmon=month+1
  if month gt 0 then begin
    tdy=tdy-total(months_days[0:month-1])
  endif
  
  RETURN, [tyr,tmon,tdy]
END