program main use datatype_module implicit none include 'netcdf.inc' include '/usr/include/hdf/dffunc.inc' c define variables integer*4 latmax,lonmax,daymax,binmax,itmax integer*4 ilat,ilon integer*4 ical,imod,ipar real :: pi integer it,prof,ibin,ibin1,imon,imon1 integer numfichpar,numfichmodis,numfichmod,numfichpara integer ifile,ifile3 integer ifile4,ifile5,ifile7 integer ncid,varid,dimid,status,rcode,dir integer tart(3),count(3),vdims(3) integer nctype,nvdim,nvatts,lenstr,j,ndsize character*1024 strbuf character*50 FILE_NAME integer retval integer lat_varid,lon_varid integer refcloudmod2_varid,refcloudmod1_varid,refcloudpar_varid integer lat_dimid, lon_dimid integer dimids(2) integer ndims parameter (ndims=2) parameter (latmax=90) parameter (lonmax=180) parameter (daymax=31) parameter (binmax=502) parameter (itmax=47304000) !365*2*lonmax*latmax real*4 lonmod(lonmax),latmod(latmax) real*4 bin(binmax) character*125 file3,file5,file4,file7 integer*2 jour character*2 day character*2 month character montha character*2 monthb integer*2 mois,moisa real*8,allocatable :: temps(:) real*4,allocatable :: lat(:) real*4,allocatable :: lon(:) integer*1,allocatable :: clf(:) integer*2,allocatable :: cf(:) real*4,allocatable :: cft(:) integer*2,allocatable :: refl(:) real*4,allocatable :: ref(:) integer*1,allocatable :: fisp(:) integer*2,allocatable :: fis(:) integer*2,allocatable :: thetasp(:,:) integer*2,allocatable :: thetavp(:,:) real*4,allocatable :: thetas(:,:) real*4,allocatable :: thetav(:,:) integer*2,allocatable :: fip(:,:) real*4,allocatable :: fi(:,:) real*4,allocatable :: fiv(:,:) integer*1,allocatable :: cosa(:,:) integer*2,allocatable :: rad865p(:,:) real*4,allocatable :: rad865(:,:) real*4,allocatable :: scat_ang(:,:) real*4,allocatable :: rad(:),theta(:),refle(:) real*8 reflat,reflon real*8 modlat,modlon integer*4 sd_id integer*4 dims integer*4,allocatable :: start(:),stride(:) integer*4 sds_id,sds_index character*64 sds_name character*256 dimlist integer*4 rank,datatype integer*4 dimv(17) real*4 ::cltcalipso(lonmax,latmax) real*4 ::cfcal(latmax,lonmax) real*4 :: ind_cal(latmax,lonmax),refmod_cal(itmax), &refpar_cal(itmax), refmod_mod(itmax),refpar_mod(itmax), &refmod_par(itmax),refpar_par(itmax),refmod250_mod(itmax), &refmod250_par(itmax),refmod250_cal(itmax) integer*4 iilat,iilon real*8 parlat,parlon,x,y real*4 refpar(latmax,lonmax),indrpar(latmax,lonmax), &refmod(latmax,lonmax),refcloudmod2(latmax,lonmax), &indmodiscloud(latmax,lonmax),refclearmod2(latmax,lonmax), &indmodisclear(latmax,lonmax),refcloudpar(latmax,lonmax), &indparcloud(latmax,lonmax),refcloudmod1(latmax,lonmax), &indmodis1cloud(latmax,lonmax) real*4 crefmod2_mon(latmax,lonmax),indmod2_mon(latmax,lonmax), &crefmod1_mon(latmax,lonmax),indmod1_mon(latmax,lonmax), & crefpar_mon(latmax,lonmax),indpar_mon(latmax,lonmax) real*8 indrmod(latmax,lonmax), &indrmod250(latmax,lonmax),refmod250(latmax,lonmax) real*8 refhpar(latmax,lonmax,binmax),refhmod(latmax,lonmax,binmax) & ,refhmod2(latmax,lonmax,binmax),refparcdf(latmax,lonmax,binmax), & refmodcdf(latmax,lonmax,binmax),refmod2cdf(latmax,lonmax,binmax) real*8 pourc_oce(latmax,lonmax),pour_oce(latmax,lonmax) real*8 refl_scale integer attr_index bin(1)=0. do ibin=2,binmax bin(ibin)=int(bin(ibin-1))+2. enddo do ibin=1,binmax bin(ibin)=bin(ibin)/1000. ! write(*,*) 'bin', bin(ibin) enddo ical=0. imod=0. ipar=0. pi=acos(-1.0) open(13,file= &'/homedata/dkonsta/modele/grilles_lmdz/longitude_2x2') do ilon=1,lonmax read(13,*) lonmod(ilon) !lmdz enddo open(14,file= &'/homedata/dkonsta/modele/grilles_lmdz/latitude_2x2') do ilat=1,latmax read(14,*) latmod(ilat) !lmdz latmod(ilat)=-latmod(ilat) enddo open(50,file= & '/homedata/dkonsta/modele/grilles_lmdz/land_ocean_mask2_2d') do ilat=1,latmax do ilon=1,lonmax read(50,*) pour_oce(ilat,ilon) enddo enddo do ilat=1,latmax do ilon=1,lonmax pourc_oce(ilat,ilon)=pour_oce(latmax+1-ilat,ilon) enddo enddo numfichpar=0. numfichpara=0. numfichmod=0. numfichmodis=0. ! CF CALIPSO open(5,file= & '/homedata/dkonsta/lidar/listfichiers/list200708_day') 885 read(5,100,end=995) file5 month=file5(76:77) read(month,'i2')mois do ilon=1,lonmax do ilat=1,latmax cfcal(ilat,ilon)=0. refpar(ilat,ilon)=0. refmod(ilat,ilon)=0. refmod250(ilat,ilon)=0. indrpar(ilat,ilon)=0. indrmod(ilat,ilon)=0. indrmod250(ilat,ilon)=0. refcloudmod2(ilat,ilon)=0. indmodiscloud(ilat,ilon)=0. refclearmod2(ilat,ilon)=0. indmodisclear(ilat,ilon)=0. refcloudpar(ilat,ilon)=0. indparcloud(ilat,ilon)=0. refcloudmod1(ilat,ilon)=0. indmodis1cloud(ilat,ilon)=0. do ibin=1,binmax refhpar(ilat,ilon,ibin)=0. refhmod(ilat,ilon,ibin)=0. refhmod2(ilat,ilon,ibin)=0. refparcdf(ilat,ilon,ibin)=0. refmodcdf(ilat,ilon,ibin)=0. refmod2cdf(ilat,ilon,ibin)=0. enddo enddo enddo !open calipso write(*,*) 'read file cf calipso', file5 ncid=ncopn(file5,ncnowrit,rcode) status=nf_inq_varid(ncid,'cltcalipso',varid) call ncvinq(ncid,varid,strbuf,nctype,nvdim,vdims,nvatts,rcode) lenstr=1 do j=1,nvdim call ncdinq(ncid,vdims(j),strbuf,ndsize,rcode) lenstr=lenstr*ndsize tart(j)=1 count(j)=ndsize enddo call ncvgt(ncid,varid,tart,count,cltcalipso,rcode) do ilon=1,lonmax do ilat=1,latmax if (cltcalipso(ilon,ilat).ge.0.) then cfcal(ilat,ilon)=cltcalipso(ilon,ilat) else cfcal(ilat,ilon)=-9999. endif enddo enddo status=nf_close(ncid) ! reflectance MODIS 100 format(A120) open(3,file= &'/homedata/dkonsta/modis/listfichiers/listmodis_200708_day') do ifile3=1,10000 883 read(3,100,end=993) file3 if (file3(82:83).ne.file5(74:75).or.file3(85:86) & .ne.file5(76:77).or.file3(88:89).ne.file5(78:79)) cycle write(*,*) 'read file ref modis', file3 !ouvrir modis numfichmod=numfichmod+1 if (numfichmod.gt.1) then deallocate (temps,lat,lon,refl,ref) endif c open hdf file sd_id=sfstart(file3,DFACC_READ) if (sd_id.eq.-1)then write(*,*) 'Error sfstart' stop endif c latitude !valeur faux: -Inf sds_index=0 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dims,datatype,dimlist) ! write(*,*) 'attributes',rank,dims,sds_name allocate(lat(dims)) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 ! write(*,*) 'status',sds_id,start,stride,dims status= sfrdata(sds_id,start,stride,dims,lat) status=sfendacc(sds_id) ! write(*,*) 'lat=', lat(150:190) deallocate (start,stride) c longitude !valeur faux: -Inf sds_index=1 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dims,datatype,dimlist) allocate(lon(dims)) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 status=sfrdata(sds_id,start,stride,dims,lon) status=sfendacc(sds_id) deallocate (start,stride) c time !valeur faux : -Inf sds_index=2 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dims,datatype,dimlist) ! write(*,*) 'attributes',rank,dims,sds_name allocate(temps(dims)) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 ! write(*,*) 'status',sds_id,start,stride,dims status=sfrdata(sds_id,start,stride,dims,temps) ! write(*,*) 'temps=', temps(1:10) status=sfendacc(sds_id) deallocate (start,stride) c refl !valeur faux : 65535 sds_index=6 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dims,datatype,dimlist) if (sds_name.eq.'MYD021KM_EV_250_Aggr1km_RefSB_Band2') then ! write(*,*) 'attributes',rank,dims,sds_name allocate(refl(dims),ref(dims)) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 ! write(*,*) 'status',sds_id,start,stride,dims status=sfrdata(sds_id,start,stride,dims,refl) ! write(*,*) 'refl=', refl(1:10) status=sfendacc(sds_id) deallocate (start,stride) c reflectance scale !valeur faux : 65535 sds_index=6 sds_id=sfselect(sd_id,sds_index) attr_index=sffattr(sds_id,'reflectance_scale') status=sfrnatt(sds_id,attr_index,refl_scale) ! write(*,*) 'scale',refl_scale,attr_index endif prof=dims status = sfend (sd_id); if (sds_name.ne.'MYD021KM_EV_250_Aggr1km_RefSB_Band2') then allocate(ref(prof),refl(prof)) do it=1,prof refl(it)=-9999. enddo endif do it=1,prof ref(it)=refl(it) if(ref(it).gt.0.) then ref(it)=ref(it)*refl_scale else ref(it)=-9999. endif reflat=1000 reflon=1000 do iilat=1,latmax ! cherche la boite de lat modelisee correspondante x=abs(lat(it)-latmod(iilat)) if (x.lt.reflat) then ilat=iilat reflat=x endif enddo c do iilon=1,lonmax ! cherche boite lon modelisée y=abs(lon(it)-lonmod(iilon)) if (y.lt.reflon) then ilon=iilon reflon=y endif enddo if(ref(it).ne.-9999.) then refmod(ilat,ilon)=refmod(ilat,ilon)+ref(it) indrmod(ilat,ilon)=indrmod(ilat,ilon)+1 do ibin=1,binmax-1 if (bin(ibin).le.ref(it). & and.ref(it).lt.bin(ibin+1))then refhmod(ilat,ilon,ibin)=refhmod(ilat,ilon,ibin)+1. endif enddo endif enddo ! REFLECTANCE PARASOL open(4,file= &'/homedata/dkonsta/parasol/listfichiers/listparascoloc_200708') do ifile4=1,10000 884 read(4,100,end=994) file4 if (file3(80:98).ne.file4(78:96)) cycle write(*,*) 'read file ref parasol', file4 !ouvrir parasol numfichpara=numfichpara+1 if (numfichpara.gt.1.) then deallocate (fis,fisp,fiv,thetas,fi,thetav, & cosa,rad865,rad865p,rad,scat_ang,thetasp,thetavp,fip,theta,refle) endif deallocate(temps,lat,lon) c open hdf file sd_id=sfstart(file4,DFACC_READ) if (sd_id.eq.-1)then write(*,*) 'Error sfstart' stop endif c latitude !valeur faux: -Inf sds_index=0 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dimv,datatype,dimlist) ! write(*,*) 'attributes',rank,dimv,sds_name allocate(lat(dimv(1))) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 ! write(*,*) 'status',sds_id,start,stride,dimv status= sfrdata(sds_id,start,stride,dimv,lat) ! write(*,*) 'lat=', lat(1150:1190) status=sfendacc(sds_id) deallocate (start,stride) c longitude !valeur faux: -Inf sds_index=1 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dimv,datatype,dimlist) allocate(lon(dimv(1))) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 status=sfrdata(sds_id,start,stride,dimv,lon) status=sfendacc(sds_id) deallocate (start,stride) c time !valeur faux : -Inf sds_index=2 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dimv,datatype,dimlist) ! write(*,*) 'attributes',rank,dimv,sds_name allocate(temps(dimv(1))) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 ! write(*,*) 'status',sds_id,start,stride,dimv status=sfrdata(sds_id,start,stride,dimv,temps) ! write(*,*) 'temps=', temps(1:10) status=sfendacc(sds_id) deallocate (start,stride) c solar_azimuth_angle !valeur faux: 255 (probleme) sds_index=5 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dimv,datatype,dimlist) ! write(*,*) 'attributes',rank,dimv,sds_name allocate(fisp(dimv(1)),fis(dimv(1))) allocate(start(rank),stride(rank)) start(:)=0. stride(:)=1. ! write(*,*) 'status',sds_id,start,stride,dimv status=sfrdata(sds_id,start,stride,dimv,fisp) ! write(*,*) 'fisp=', fisp(1:10) status=sfendacc(sds_id) deallocate (start,stride) c solar zenith angle !valeur faux :65535 sds_index=7 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dimv,datatype,dimlist) ! write(*,*) 'attributes',rank,dimv,sds_name allocate(thetas(dimv(1),dimv(2))) allocate(thetasp(dimv(1),dimv(2))) allocate(scat_ang(dimv(1),dimv(2))) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 ! write(*,*) 'status',sds_id,start,stride,dimv status= sfrdata(sds_id,start,stride,dimv(1:rank),thetasp) ! write(*,*) 'solarzenithangle=', (thetas(10000,j),j=1,dimv(1)) status=sfendacc(sds_id) deallocate (start,stride) prof=dimv(2) dir=dimv(1) c view zenith angle ! valeur faux: 65535 sds_index=8 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dimv,datatype,dimlist) ! write(*,*) 'attributes',rank,dimv(1),dimv(2),sds_name allocate(thetav(dimv(1),dimv(2))) allocate(thetavp(dimv(1),dimv(2))) allocate(theta(dimv(2))) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 ! write(*,*) 'status',sds_id,start,stride status= sfrdata(sds_id,start,stride,dimv(1:rank),thetavp) ! write(*,*) 'view_angle=', (thetav(10000,j),j=1,dimv(1)) status=sfendacc(sds_id) deallocate (start,stride) c relative azimuth angle ! valeur faux:65535 sds_index=9 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dimv,datatype,dimlist) ! write(*,*) 'attributes',rank,dimv(1),dimv(2),sds_name allocate(fi(dimv(1),dimv(2)),fiv(dimv(1),dimv(2))) allocate(fip(dimv(1),dimv(2))) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 ! write(*,*) 'status',sds_id,start,stride status= sfrdata(sds_id,start,stride,dimv(1:rank),fip) ! write(*,*) 'ralative azimuth=', (fi(10000,j),j=1,dimv(1)) status=sfendacc(sds_id) deallocate (start,stride) c rdelta thetav cos !valeur faux: -128 sds_index=10 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dimv,datatype,dimlist) ! write(*,*) 'attributes',rank,dimv(1),dimv(2),sds_name allocate(cosa(dimv(1),dimv(2))) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 ! write(*,*) 'status',sds_id,start,stride status= sfrdata(sds_id,start,stride,dimv(1:rank),cosa) ! write(*,*) 'ralative azimuth=', (cosa(10000,j),j=1,dimv(1)) status=sfendacc(sds_id) deallocate (start,stride) c radiance 865 ! valeur faux: -32768 sds_index=19 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dimv,datatype,dimlist) ! write(*,*) 'attributes',rank,dimv(1),dimv(2),sds_name allocate(rad865p(dimv(1),dimv(2)),rad(dimv(2)),refle(dimv(2))) allocate(rad865(dimv(1),dimv(2))) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 ! write(*,*) 'status',sds_id,start,stride status= sfrdata(sds_id,start,stride,dimv(1:rank),rad865p) ! write(*,*) 'rad865=', (rad865(10000,j),j=1,dimv(1)) status=sfendacc(sds_id) deallocate (start,stride) status = sfend (sd_id); do it=1,prof refle(it)=-9999. fis(it)=fisp(it) if (fis(it).lt.0.) then fis(it)=fis(it) + 256. endif fis(it)=fis(it)*1.4 do j=1,dir fi(j,it)=fip(j,it) if (fi(j,it).lt.0.) then fi(j,it)=fi(j,it)+65536. endif rad865(j,it)=rad865p(j,it) thetav(j,it)=thetavp(j,it) thetas(j,it)=thetasp(j,it) if (thetav(j,it).lt.0.) then thetav(j,it)=thetav(j,it)+65536. endif if (thetas(j,it).lt.0.) then ! write(*,*) 'thetaprin', thetas(j,it) thetas(j,it)=thetas(j,it)+65536. ! write(*,*) 'thetameta', thetas(j,it) endif thetav(j,it)=thetav(j,it)*0.0015 thetas(j,it)=thetas(j,it)*0.0015 fi(j,it)=fi(j,it)*0.006 cosa(j,it)=cosa(j,it)*0.0016 rad865(j,it)=rad865(j,it)*0.0001 fiv(j,it)=fis(it)-fi(j,it) if (thetas(j,it).eq.98.30250.or.thetav(j,it).eq.98.30250.or. & fis(it).eq.357.or.fi(j,it).eq.393.21) then scat_ang(j,it)=-9999. else scat_ang(j,it)=acos(cos(thetas(j,it)*pi/180.)*cos(thetav(j,it)* & pi/180.)+sqrt(1-cos(thetas(j,it)*pi/180.)* & cos(thetas(j,it)*pi/180.))*sqrt(1-cos(thetav(j,it)*pi/180.)* & cos(thetav(j,it)*pi/180.))*cos(-fi(j,it)*pi/180.))*180./pi endif if ((thetav(j,it).gt.24.5.and.thetav(j,it).lt.29.5).and. & (((fi(j,it).gt.335.).and.(fi(j,it).lt.345.)).or. & ((fi(j,it).gt.15.).and.(fi(j,it).lt.25.))).and. & (rad865(j,it).ne.-3.2767)) then rad(it)=rad865(j,it) theta(it)=thetas(j,it) refle(it)=rad(it)/cos(theta(it)*pi/180.) endif enddo reflat=1000 reflon=1000 do iilat=1,latmax ! cherche la boite de lat modelisee correspondante x=abs(lat(it)-latmod(iilat)) if (x.lt.reflat) then ilat=iilat reflat=x endif enddo c do iilon=1,lonmax ! cherche boite lon modelisée y=abs(lon(it)-lonmod(iilon)) if (y.lt.reflon) then ilon=iilon reflon=y endif enddo if(refle(it).ge.0.) then refpar(ilat,ilon)=refpar(ilat,ilon)+refle(it) indrpar(ilat,ilon)=indrpar(ilat,ilon)+1 do ibin=1,binmax-1 if (bin(ibin).le.refle(it). & and.refle(it).lt.bin(ibin+1))then refhpar(ilat,ilon,ibin)=refhpar(ilat,ilon,ibin)+1. endif enddo endif enddo ! reflectance MODIS 250m open(7,file= &'/homedata/dkonsta/modis/listfichiers/listmodis250m_200708_day') do ifile7=1,10000 887 read(7,100,end=997) file7 if (file3(80:98).ne.file7(87:105)) cycle write(*,*) 'read file ref modis250m', file7 !ouvrir modis deallocate (temps,lat,lon,refl,ref) c open hdf file sd_id=sfstart(file7,DFACC_READ) if (sd_id.eq.-1)then write(*,*) 'Error sfstart' stop endif c latitude !valeur faux: -Inf sds_index=0 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dims,datatype,dimlist) ! write(*,*) 'attributes',rank,dims,sds_name allocate(lat(dims)) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 ! write(*,*) 'status',sds_id,start,stride,dims status= sfrdata(sds_id,start,stride,dims,lat) status=sfendacc(sds_id) ! write(*,*) 'lat=', lat(150:190) deallocate (start,stride) c longitude !valeur faux: -Inf sds_index=1 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dims,datatype,dimlist) allocate(lon(dims)) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 status=sfrdata(sds_id,start,stride,dims,lon) status=sfendacc(sds_id) deallocate (start,stride) c time !valeur faux : -Inf sds_index=2 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dims,datatype,dimlist) ! write(*,*) 'attributes',rank,dims,sds_name allocate(temps(dims)) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 ! write(*,*) 'status',sds_id,start,stride,dims status=sfrdata(sds_id,start,stride,dims,temps) ! write(*,*) 'temps=', temps(1:10) status=sfendacc(sds_id) deallocate (start,stride) c refl !valeur faux : 65535 sds_index=8 sds_id=sfselect(sd_id,sds_index) status=sfginfo(sds_id,sds_name,rank,dims,datatype,dimlist) ! write(*,*) 'attributes',rank,dims,sds_name,sds_id,dimlist allocate(refl(dims),ref(dims)) allocate(start(rank),stride(rank)) start(:)=0 stride(:)=1 ! write(*,*) 'status',sds_id,start,stride,dims status=sfrdata(sds_id,start,stride,dims,refl) ! write(*,*) 'refl=', refl(1:10) status=sfendacc(sds_id) deallocate (start,stride) prof=dims c reflectance scale !valeur faux : 65535 sds_index=8 sds_id=sfselect(sd_id,sds_index) attr_index=sffattr(sds_id,'reflectance_scale') status=sfrnatt(sds_id,attr_index,refl_scale) ! write(*,*) 'scale',refl_scale,attr_index status = sfend (sd_id); do it=1,prof ref(it)=refl(it) if(ref(it).gt.0.) then ref(it)=ref(it)*refl_scale else ref(it)=-9999. endif reflat=1000 reflon=1000 do iilat=1,latmax ! cherche la boite de lat modelisee correspondante x=abs(lat(it)-latmod(iilat)) if (x.lt.reflat) then ilat=iilat reflat=x endif enddo c do iilon=1,lonmax ! cherche boite lon modelisée y=abs(lon(it)-lonmod(iilon)) if (y.lt.reflon) then ilon=iilon reflon=y endif enddo if(ref(it).ne.-9999.) then refmod250(ilat,ilon)=refmod250(ilat,ilon)+ref(it) indrmod250(ilat,ilon)=indrmod250(ilat,ilon)+1 do ibin=1,binmax-1 if (bin(ibin).le.ref(it). & and.ref(it).lt.bin(ibin+1))then refhmod2(ilat,ilon,ibin)=refhmod2(ilat,ilon,ibin)+1. endif enddo endif enddo go to 887 !next file modis reflectance 250m enddo 997 continue !end of files modis reflectance 250m close(7) go to 884 !next file parasol reflectance enddo 994 continue !end of files parasol reflectance close(4) go to 883 !next file modis reflectance enddo 993 continue !end of files modis reflectance close(3) do ilon=1,lonmax do ilat=1,latmax do ibin=1,binmax if (ibin.eq.1.) then refparcdf(ilat,ilon,ibin)=refhpar(ilat,ilon,ibin) refmodcdf(ilat,ilon,ibin)=refhmod(ilat,ilon,ibin) refmod2cdf(ilat,ilon,ibin)=refhmod2(ilat,ilon,ibin) else refparcdf(ilat,ilon,ibin)= & refparcdf(ilat,ilon,ibin-1)+refhpar(ilat,ilon,ibin) refmodcdf(ilat,ilon,ibin)= & refmodcdf(ilat,ilon,ibin-1)+refhmod(ilat,ilon,ibin) refmod2cdf(ilat,ilon,ibin)= & refmod2cdf(ilat,ilon,ibin-1)+refhmod2(ilat,ilon,ibin) endif enddo enddo enddo do ilon=1,lonmax do ilat=1,latmax do ibin=1,binmax-1 if (indrpar(ilat,ilon).gt.0.) then refparcdf(ilat,ilon,ibin)= & refparcdf(ilat,ilon,ibin)/indrpar(ilat,ilon) else refparcdf(ilat,ilon,ibin)=-9999. endif if (indrmod(ilat,ilon).gt.0.) then refmodcdf(ilat,ilon,ibin)= & refmodcdf(ilat,ilon,ibin)/indrmod(ilat,ilon) else refmodcdf(ilat,ilon,ibin)=-9999. endif if (indrmod250(ilat,ilon).gt.0.) then refmod2cdf(ilat,ilon,ibin)= & refmod2cdf(ilat,ilon,ibin)/indrmod250(ilat,ilon) else refmod2cdf(ilat,ilon,ibin)=-9999. endif enddo enddo enddo do ilon=1,lonmax do ilat=1,latmax if (indrpar(ilat,ilon).gt.0.) then refpar(ilat,ilon)=refpar(ilat,ilon)/indrpar(ilat,ilon) else refpar(ilat,ilon)=-9999. endif if (indrmod(ilat,ilon).gt.0.) then refmod(ilat,ilon)=refmod(ilat,ilon)/indrmod(ilat,ilon) else refmod(ilat,ilon)=-9999. endif if (indrmod250(ilat,ilon).gt.0.) then refmod250(ilat,ilon)= & refmod250(ilat,ilon)/indrmod250(ilat,ilon) else refmod250(ilat,ilon)=-9999. endif do ibin=1,binmax if (refhmod2(ilat,ilon,ibin).ge.0.and. & refmod2cdf(ilat,ilon,ibin).ge.0.and. & (1-cfcal(ilat,ilon)).lt.refmod2cdf(ilat,ilon,ibin))then refcloudmod2(ilat,ilon)=refcloudmod2(ilat,ilon)+ & refhmod2(ilat,ilon,ibin)*(bin(ibin)+0.001) indmodiscloud(ilat,ilon)= & indmodiscloud(ilat,ilon)+refhmod2(ilat,ilon,ibin) endif if (refhpar(ilat,ilon,ibin).ge.0.and. & refparcdf(ilat,ilon,ibin).ge.0.and. & (1-cfcal(ilat,ilon)).lt.refparcdf(ilat,ilon,ibin))then refcloudpar(ilat,ilon)=refcloudpar(ilat,ilon)+ & refhpar(ilat,ilon,ibin)*(bin(ibin)+0.001) indparcloud(ilat,ilon)= & indparcloud(ilat,ilon)+refhpar(ilat,ilon,ibin) endif if (refhmod2(ilat,ilon,ibin).ge.0.and. & refmod2cdf(ilat,ilon,ibin).ge.0.and. & (1-cfcal(ilat,ilon)).gt.refmod2cdf(ilat,ilon,ibin))then refclearmod2(ilat,ilon)=refclearmod2(ilat,ilon)+ & refhmod2(ilat,ilon,ibin)*(bin(ibin)+0.001) indmodisclear(ilat,ilon)= & indmodisclear(ilat,ilon)+refhmod2(ilat,ilon,ibin) endif if (refhmod(ilat,ilon,ibin).ge.0.and. & refmodcdf(ilat,ilon,ibin).ge.0.and. & (1-cfcal(ilat,ilon)).lt.refmodcdf(ilat,ilon,ibin))then refcloudmod1(ilat,ilon)=refcloudmod1(ilat,ilon)+ & refhmod(ilat,ilon,ibin)*(bin(ibin)+0.001) indmodis1cloud(ilat,ilon)= & indmodis1cloud(ilat,ilon)+refhmod(ilat,ilon,ibin) endif enddo enddo enddo do ilon=1,lonmax do ilat=1,latmax if (indmodiscloud(ilat,ilon).gt.0.) then refcloudmod2(ilat,ilon)= & refcloudmod2(ilat,ilon)/indmodiscloud(ilat,ilon) else refcloudmod2(ilat,ilon)=-9999. endif if (indparcloud(ilat,ilon).gt.0.) then refcloudpar(ilat,ilon)= & refcloudpar(ilat,ilon)/indparcloud(ilat,ilon) else refcloudpar(ilat,ilon)=-9999. endif if (indmodis1cloud(ilat,ilon).gt.0.) then refcloudmod1(ilat,ilon)= & refcloudmod1(ilat,ilon)/indmodis1cloud(ilat,ilon) else refcloudmod1(ilat,ilon)=-9999. endif if (pourc_oce(ilat,ilon).gt.0.5) then refcloudmod2(ilat,ilon)=-9999. refcloudmod1(ilat,ilon)=-9999. refcloudpar(ilat,ilon)=-9999. endif enddo enddo c create the file FILE_NAME='CRef_'//file5(72:79)//'_CFMIP.nc' retval=nf_create('/bdd/CFMIP/GOCCP/Dimitra/CRef/daily/'//FILE_NAME & ,nf_clobber,ncid) if (retval.ne.nf_noerr)then write(*,*) 'Error sfstart' stop endif retval=nf_put_att_text(ncid, nf_global, 'Description', 17, & 'Cloud Reflectance') retval=nf_put_att_text(ncid, nf_global, 'Date', 8, & FILE_NAME(6:13)) retval=nf_put_att_text(ncid, nf_global, 'Author', 51, & 'Dimitra Konsta, Helene Chepfer, Jean-Louis Dufresne') retval=nf_put_att_text(ncid, nf_global, 'references', 55, & 'http://climserv.ipsl.polytechnique.fr/cfmip-atrain.html') c define dimensions retval = nf_def_dim(ncid, 'latitude', latmax, lat_dimid) retval = nf_def_dim(ncid, 'longitude', lonmax, lon_dimid) c define the coordinate variables retval = nf_def_var(ncid, 'latitude', NF_REAL, 1, lat_dimid, & lat_varid) retval = nf_def_var(ncid, 'longitude', NF_REAL, 1, lon_dimid, & lon_varid) c assign units attributes to var data retval = nf_put_att_text(ncid, lat_varid, 'lat_name', 8, & 'latitude') retval = nf_put_att_text(ncid, lat_varid, 'units', 13, & 'degrees_north') retval = nf_put_att_text(ncid, lon_varid, 'lon_name', 9, & 'longitude') retval = nf_put_att_text(ncid, lon_varid, 'units', 12, & 'degrees_east') c define the netcdf variables dimids(1)= lat_dimid dimids(2)= lon_dimid retval = nf_def_var(ncid, 'CRef_mod250', NF_REAL, ndims, dimids, & refcloudmod2_varid) retval = nf_def_var(ncid, 'CRef_mod1km', NF_REAL, ndims, dimids, & refcloudmod1_varid) retval = nf_def_var(ncid, 'CRef_par', NF_REAL, ndims, dimids, & refcloudpar_varid) c assign untis attributes to netcdf variables retval = nf_put_att_text(ncid, refcloudmod2_varid, 'name', 28, & 'Cloud reflectance MODIS 250m') retval = nf_put_att_text(ncid, refcloudmod2_varid, 'units',14, & 'Arbitrary unit') retval = nf_put_att_text(ncid,refcloudmod2_varid,'missing value', & 5, '-9999') retval = nf_put_att_text(ncid, refcloudmod1_varid, 'name', 27, & 'Cloud reflectance MODIS 1km') retval = nf_put_att_text(ncid, refcloudmod1_varid, 'units',14, & 'Arbitrary unit') retval = nf_put_att_text(ncid,refcloudmod1_varid,'missing value', & 5, '-9999') retval = nf_put_att_text(ncid, refcloudpar_varid, 'name', 25, & 'Cloud reflectance PARASOL') retval = nf_put_att_text(ncid, refcloudpar_varid, 'units',14, & 'Arbitrary unit') retval = nf_put_att_text(ncid,refcloudpar_varid,'missing value', & 5, '-9999') c end define mode retval = nf_enddef(ncid) c write the viariable data retval = nf_put_var_real(ncid, lat_varid, latmod) if (retval.ne.nf_noerr)then write(*,*) 'nf_put_var_real' stop endif retval = nf_put_var_real(ncid, lon_varid, lonmod) retval = nf_put_var_real(ncid, refcloudmod2_varid, refcloudmod2) retval = nf_put_var_real(ncid, refcloudmod1_varid, refcloudmod1) retval = nf_put_var_real(ncid, refcloudpar_varid, refcloudpar) c close the file retval = nf_close(ncid) if (retval.ne.nf_noerr)then write(*,*) 'nf_close' stop endif if(mois.eq.imon1) then imon=imon do ilon=1,lonmax do ilat=1,latmax if (refcloudmod2(ilat,ilon).ge.0.) then crefmod2_mon(ilat,ilon)= & crefmod2_mon(ilat,ilon)+refcloudmod2(ilat,ilon) indmod2_mon(ilat,ilon)=indmod2_mon(ilat,ilon)+1. endif if (refcloudmod1(ilat,ilon).ge.0.) then crefmod1_mon(ilat,ilon)= & crefmod1_mon(ilat,ilon)+refcloudmod1(ilat,ilon) indmod1_mon(ilat,ilon)=indmod1_mon(ilat,ilon)+1. endif if (refcloudpar(ilat,ilon).ge.0.) then crefpar_mon(ilat,ilon)= & crefpar_mon(ilat,ilon)+refcloudpar(ilat,ilon) indpar_mon(ilat,ilon)=indpar_mon(ilat,ilon)+1. endif enddo enddo else do ilon=1,lonmax do ilat=1,latmax imon=imon+1. if (indmod2_mon(ilat,ilon).gt.0) then crefmod2_mon(ilat,ilon)= & crefmod2_mon(ilat,ilon)/indmod2_mon(ilat,ilon) else crefmod2_mon(ilat,ilon)=-9999. endif if (indmod1_mon(ilat,ilon).gt.0) then crefmod1_mon(ilat,ilon)= & crefmod1_mon(ilat,ilon)/indmod1_mon(ilat,ilon) else crefmod1_mon(ilat,ilon)=-9999. endif if (indpar_mon(ilat,ilon).gt.0) then crefpar_mon(ilat,ilon)= & crefpar_mon(ilat,ilon)/indpar_mon(ilat,ilon) else crefpar_mon(ilat,ilon)=-9999. endif enddo enddo moisa= mois-1. c create the file if (moisa.lt.10.) then write(montha,'i1') moisa FILE_NAME='CRef_'//file5(72:75)//'0'//montha//'_CFMIP.nc' else write(monthb,'i2') moisa FILE_NAME='CRef_'//file5(72:75)//monthb//'_CFMIP.nc' endif retval=nf_create('/bdd/CFMIP/GOCCP/Dimitra/CRef/'//FILE_NAME & ,nf_clobber,ncid) if (retval.ne.nf_noerr)then write(*,*) 'Error sfstart' stop endif retval=nf_put_att_text(ncid, nf_global, 'Description', 17, & 'Cloud Reflectance') retval=nf_put_att_text(ncid, nf_global, 'Date', 8, & FILE_NAME(6:13)) retval=nf_put_att_text(ncid, nf_global, 'Author', 51, & 'Dimitra Konsta, Helene Chepfer, Jean-Louis Dufresne') retval=nf_put_att_text(ncid, nf_global, 'references', 55, & 'http://climserv.ipsl.polytechnique.fr/cfmip-atrain.html') c define dimensions retval = nf_def_dim(ncid, 'latitude', latmax, lat_dimid) retval = nf_def_dim(ncid, 'longitude', lonmax, lon_dimid) c define the coordinate variables retval = nf_def_var(ncid, 'latitude', NF_REAL, 1, lat_dimid, & lat_varid) retval = nf_def_var(ncid, 'longitude', NF_REAL, 1, lon_dimid, & lon_varid) c assign units attributes to var data retval = nf_put_att_text(ncid, lat_varid, 'lat_name', 8, & 'latitude') retval = nf_put_att_text(ncid, lat_varid, 'units', 13, & 'degrees_north') retval = nf_put_att_text(ncid, lon_varid, 'lon_name', 9, & 'longitude') retval = nf_put_att_text(ncid, lon_varid, 'units', 12, & 'degrees_east') c define the netcdf variables dimids(1)= lat_dimid dimids(2)= lon_dimid retval = nf_def_var(ncid, 'CRef_mod250', NF_REAL, ndims, dimids, & refcloudmod2_varid) retval = nf_def_var(ncid, 'CRef_mod1km', NF_REAL, ndims, dimids, & refcloudmod1_varid) retval = nf_def_var(ncid, 'CRef_par', NF_REAL, ndims, dimids, & refcloudpar_varid) c assign untis attributes to netcdf variables retval = nf_put_att_text(ncid, refcloudmod2_varid, 'name', 28, & 'Cloud reflectance MODIS 250m') retval = nf_put_att_text(ncid, refcloudmod2_varid, 'units',14, & 'Arbitrary unit') retval = nf_put_att_text(ncid,refcloudmod2_varid,'missing value', & 5, '-9999') retval = nf_put_att_text(ncid, refcloudmod1_varid, 'name', 27, & 'Cloud reflectance MODIS 1km') retval = nf_put_att_text(ncid, refcloudmod1_varid, 'units',14, & 'Arbitrary unit') retval = nf_put_att_text(ncid,refcloudmod1_varid,'missing value', & 5, '-9999') retval = nf_put_att_text(ncid, refcloudpar_varid, 'name', 25, & 'Cloud reflectance PARASOL') retval = nf_put_att_text(ncid, refcloudpar_varid, 'units',14, & 'Arbitrary unit') retval = nf_put_att_text(ncid,refcloudpar_varid,'missing value', & 5, '-9999') c end define mode retval = nf_enddef(ncid) c write the viariable data retval = nf_put_var_real(ncid, lat_varid, latmod) if (retval.ne.nf_noerr)then write(*,*) 'nf_put_var_real' stop endif retval = nf_put_var_real(ncid, lon_varid, lonmod) retval = nf_put_var_real(ncid, refcloudmod2_varid, crefmod2_mon) retval = nf_put_var_real(ncid, refcloudmod1_varid, crefmod1_mon) retval = nf_put_var_real(ncid, refcloudpar_varid, crefpar_mon) c close the file retval = nf_close(ncid) if (retval.ne.nf_noerr)then write(*,*) 'nf_close' stop endif do ilon=1,lonmax do ilat=1,latmax if (refcloudmod2(ilat,ilon).gt.-100.) then crefmod2_mon(ilat,ilon)=refcloudmod2(ilat,ilon) indmod2_mon(ilat,ilon)=1. else crefmod2_mon(ilat,ilon)=0. indmod2_mon(ilat,ilon)=0. endif if (refcloudmod1(ilat,ilon).gt.-100.) then crefmod1_mon(ilat,ilon)=refcloudmod1(ilat,ilon) indmod1_mon(ilat,ilon)=1. else crefmod1_mon(ilat,ilon)=0. indmod1_mon(ilat,ilon)=0. endif if (refcloudpar(ilat,ilon).gt.-100.) then crefpar_mon(ilat,ilon)=refcloudpar(ilat,ilon) indpar_mon(ilat,ilon)=1. else crefpar_mon(ilat,ilon)=0. indpar_mon(ilat,ilon)=0. endif enddo enddo endif imon1=mois go to 885 !next calipso file 995 continue close(5) end