subroutine inicdf !----------------------------------------------------------------------------------- ! Initializes the netcdf input files that remain open during all the chimere run ! Longitude and latitude must be read BEFORE inigeom. !----------------------------------------------------------------------------------- use netcdf use chimere_common use gridvars, only : xlati_entire_domain, xlong_entire_domain use ncvar implicit none #define NCERR(lnum) if(ncstat/=NF90_NOERR) call nc_err(ncstat,lnum,__FILE__) !----------------------------------------------------------------------------------- ! local scalars integer :: ncstat, titleLength, k, dl,we,sn,ifn_vcoo,ive,ierr integer :: tmp, nsoils_dimid, land_cat_dimid, lusef_vid ! for file exdomout.nc integer :: metro_time_dimid, metro_date_dimid, metro_zonal_dimid, metro_merid_dimid integer :: metro_layers_dimid, metro_times integer :: idstart, idyst, ihstart, imonthstart, iystart, nh ! functions integer :: idaytype ! Netcdf id integer :: metro_ncid character (len = 200) :: title, descriptlu type(ncfile) :: nc_geog ! geog file structures containing the coordinates lon, lat !---------------------------------------------------------- print*,achar(27)//'[31m o [inicdf] '//achar(27)//'[0m Open input NetCDF files' !---------------------------------------------------------- if (chimverb.gt.2) print*,' Read domain file: '//trim(fnigeo) ! Read longitude/latitude. Coordinates are read from geog file. call nc_get_struct(fnigeo, nc_geog) ! Use LANDUSEF to know if we are using USGS or IGBP-MODIS landuse ncstat=nf90_inq_varid(nc_geog%id,'LANDUSEF',lusef_vid) ncstat=nf90_inquire_attribute(nc_geog%id,lusef_vid,"description",len=titleLength) if (ncstat == nf90_noerr .and. len(title) >= titleLength) then ncstat = nf90_get_att(nc_geog%id, lusef_vid, "description", descriptlu) if (ncstat == nf90_noerr) then do k=1,190 if(descriptlu(k:k+3)=='USGS')then iludatabase=1 exit endif if(descriptlu(k:k+4)=='MODIS')then iludatabase=2 exit endif enddo endif endif print*,' Landuse characteristics found: ',trim(descriptlu),' with iludatabase=',iludatabase ! Read output grid call nc_get_coord(nc_geog) xlong_entire_domain=nc_geog%lon xlati_entire_domain=nc_geog%lat ! Looking for indices for mineral dust emissions ncstat=nf90_inq_dimid(nc_geog%id, 'soil_cat', nsoils_dimid) ncstat=nf90_inquire_dimension(nc_geog%id, nsoils_dimid, len = nsoils) ncstat=nf90_inq_dimid(nc_geog%id, 'land_cat', land_cat_dimid) ncstat=nf90_inquire_dimension(nc_geog%id, land_cat_dimid, len = land_cat) call nc_get_var(nc_geog,'SOILCTOP',SOILrelsurf_entire_domain) call nc_get_var(nc_geog,'GREENFRAC',greenfrac_entire_domain) call nc_get_var(nc_geog,'LANDUSEF',landusef_entire_domain) call nc_get_var(nc_geog,'LANDMASK',landmask_ed) ncstat=nf90_get_att(nc_geog%id,NF90_GLOBAL,'ISWATER',iswater) ncstat=nf90_get_att(nc_geog%id,NF90_GLOBAL,'ISLAKE',islake) ncstat=nf90_get_att(nc_geog%id,NF90_GLOBAL,'ISICE',isice) ncstat=nf90_get_att(nc_geog%id,NF90_GLOBAL,'ISURBAN',isurban) ncstat=nf90_get_att(nc_geog%id,NF90_GLOBAL,'ISOILWATER',isoilwater) if(chimverb.ge.2)then print*,'inicdf: geog file: iswater=',iswater print*,' islake=',islake print*,' isice=',isice print*,' isurban=',isurban print*,' isoilwater=',isoilwater endif ncstat=nf90_inq_varid(nc_geog%id,'VAR_SSO', tmp) if (ncstat.eq.NF90_NOERR) then call nc_get_var(nc_geog,'VAR_SSO',varsso_entire_domain) else if (.not.allocated(varsso_entire_domain)) & & allocate( varsso_entire_domain(nzonal_entire_domain,nmerid_entire_domain) ) varsso_entire_domain=15000. ! Mean value of maximum orographic variance endif maxvoro=maxval(varsso_entire_domain) if(chimverb.ge.2)then print*,'inicdf: successfully read and closed geog file' endif call nc_close(nc_geog) !------------------------------------------------------------------------------------ ! Initialization of time variables ! definition of 'idate' from the starting date 'idatestart' ! idate = the YYYYMMDDHH date at all times ! idtyp = the day type at all times if(.not.allocated(imonth))allocate(imonth(0:nhourrun)) if(.not.allocated(ihour)) allocate(ihour (0:nhourrun)) if(.not.allocated(idtyp)) allocate(idtyp (0:nhourrun)) if(.not.allocated(idate)) allocate(idate (0:nhourrun)) if(.not.allocated(iday)) allocate(iday (0:nhourrun)) if(.not.allocated(iyear)) allocate(iyear (0:nhourrun)) call ddate(idatestart,iystart,imonthstart,idstart,ihstart) call rdate(iystart,1,1,0,idyst) djul = interdat(idyst,idatestart)/24. do nh=0,nhourrun call reldat(idatestart,nh,idate(nh)) call ddate(idate(nh),iyear(nh),imonth(nh),iday( nh),ihour( nh)) idtyp(nh) = idaytype(idate(nh)) enddo iminute=0 isecond=0 ! For zenith ! idrefe: 1st january 2000 (for zenith based on jacobson parameterization) ! y2ktim: variable Njd in Jacobson (must be in days) call rdate(2000,01,01,00,idrefe) y2ktim = real(interdat(idrefe,idatestart))/24d0 ! Old version of zenithal angle calculation call rdate(iystart,12,21,12,idrefe) soltim = interdat(idrefe,idatestart) !----------------------------------------------------------------------------------- ! Opening EXDOMOUT.nc file if (trim(meteo_type)=='WRF') then call ini_prepwrf else print*,'EXDOMOUT file is open: ',fnmetraw(1:len_trim(fnmetraw)) ncstat=nf90_open(fnmetraw,NF90_NOWRITE,metro_ncid) ! Check correctness of dimensions ! IDs ncstat=nf90_inq_dimid(metro_ncid, 'Time', metro_time_dimid) ncstat=nf90_inq_dimid(metro_ncid, 'DateStrLen', metro_date_dimid) ncstat=nf90_inq_dimid(metro_ncid, 'west_east', metro_zonal_dimid) ncstat=nf90_inq_dimid(metro_ncid, 'south_north', metro_merid_dimid) ncstat=nf90_inq_dimid(metro_ncid, 'bottom_top', metro_layers_dimid) ! lengths ncstat=nf90_inquire_dimension(metro_ncid,metro_time_dimid,len=metro_times) if(metro_times