subroutine ini_prepwrf !------------------------------------------------------------------------------------------- ! Eric Chaxel, GIERSA, 5 March 2007 ! J.-L. Monge 2010-03-23 ! modif forecast/analysis ! modif calc_alt ! Dmitry Khvorostyanov ! 2011-09-30: continuous WRF runs treatment ! September 2016: adaptable to WRF files ! Augustin Colette May 2013 ! prepwrf: 3h WRF output fields ! new precip decumulation => no need for previous WRF file ! Sylvain Mailler, December 2014 ! New interpolation scheme for winds : spherical geometry to be able ! to use prepwrf in polar areas ! Regis Briand, 2015: prepwrf put in CHIMERE core ! This program reads netCDF wrfout file(s) from WRF, ! extracts the variables required by CHIMERE, ! performs an horizontal interpolation on the Chimere mesh, ! and generates a netCDF file. ! The program also performs the rotation of wind velocity vectors ! ! When in forecast mode, decumulation of the cumulated fields is performed assuming ! these fields are null at initial state. ! When in analysis mode, the decumulation is performed using fields from hour 23:00 ! of the previous WRF run !------------------------------------------------------------------------------------------- use netcdf use chimere_common use gridvars, only : xlati_entire_domain, xlong_entire_domain use ncvar use interp_tools implicit none #define NCERR(lnum) if(ncstat/=NF90_NOERR) call nc_err(ncstat,lnum,__FILE__) !------------------------------------------------------------------------------------------- ! Input parameters integer :: ids ! Starting date of CHIMERE simulation (YYYYMMDDHH) integer :: nhours ! Number of simulation hours namelist /args/ ids, nhours, fn_wrf, n_wrf_dates, meteo_dir character(len=*),parameter :: fniARGS = 'prepwrf.nml' ! Constant netCDF attributes character(len=*),parameter :: title='CHIMERE SUITE' character(len=*),parameter :: subtitle='Interpolated Meteo file' character(len=*),parameter :: generator='Generated by prepwrf' character(len=*),parameter :: conventions='' ! local constants character(len=*),parameter :: fn_mdom='./COORD_WRF' character(len=*),parameter :: fn_mdom_u='./COORD_WRF_U' character(len=*),parameter :: fn_mdom_v='./COORD_WRF_V' ! WRF related integer :: ifn_wrf, ifn_args ! Indices and other stuff integer :: i, nh, it, ihourcount integer :: ime,izo integer :: ierr ! array dimensions integer :: npts ! number of points in a weight file ! misc integer :: ide ! Ending date of CHIMERE simulation (YYYYMMDDHH) ! system information character(len=32) :: systime character(len=64) :: hname character(len=32) :: usrname character(len=255) :: cwd ! misc integer :: startYear, startMnth, startDay, startday1, startHour, startHour1, startMin, startSec integer :: endYear, endMnth, endDay, endHour, endMin, endSec, numFrTimes character(len=dlen) :: wrfDatestr ! NetCDF structures from the ncvar module ! WRF file structures containing the coordinates lon, lat, lev and times type(ncfile) :: nc_WRF, nc_WRFINPUT integer, allocatable :: wrf_date(:) character(len=64), allocatable :: wrf_file(:) character(len=dlen),external :: numeric2mm5date ! Functions logical :: test_wrf_grid !------------------------------------------------------------------------------------------- ! Read input arguments from namelist if (online == 0) then call opfi(ifn_args,fniARGS,'f','o') read(ifn_args,nml=args,iostat=ierr) if(ierr/=0) then call exit1(__FILE__//': problem reading arguments namelist') end if close (ifn_args) if ((fn_wrf=='').or.(fnmetraw=='').or.(ids==0).or.(nhours==0)) then call exit1('prepwrf : Not enough arguments given') end if usrname='unknown' hname='unknown' systime='unknown' cwd='unknown' ! calculate end date ide call reldat(ids,nhours,ide) if(chimverb.gt.2)then write(*,*) 'CHIMERE starting date: ids=',ids write(*,*) ' ending date: ide=',ide endif endif ! Geometric calculations for vector interpolations if (.not.allocated(chim_centers_3d))allocate(chim_centers_3d(nzonal_entire_domain,nmerid_entire_domain,3)) if (.not.allocated(chim_ulambda))allocate(chim_ulambda(nzonal_entire_domain,nmerid_entire_domain,3)) if (.not.allocated(chim_uphi))allocate(chim_uphi(nzonal_entire_domain,nmerid_entire_domain,3)) call calculate_centers(xlong_entire_domain,xlati_entire_domain, & & chim_centers_3d,chim_ulambda,chim_uphi,nzonal_entire_domain,nmerid_entire_domain) !-------------------------------- ! OFFLINE if ( online == 0 ) then if (chimverb.ge.2) print *,' Opening WRF files list '//trim(fn_wrf)//'...' call opfi(ifn_wrf, fn_wrf,'f','o') allocate( wrf_date(n_wrf_dates), wrf_file(n_wrf_dates), nh2wrf_file(0:nhourrun) ) do i = 1, n_wrf_dates read (ifn_wrf, *, iostat=ierr) wrf_date(i), wrf_file(i) enddo if (ierr/=0) then print *, '**** prepwrf : problem reading ', trim(fn_wrf) call exit1('Exiting from ini_prepwrf') endif close(ifn_wrf) ! Fill missing WRF times if any i = 1 do nh = 0, nhourrun if (i > n_wrf_dates) call exit1('*** Not enough dates in '//trim(fn_wrf)//' for date '//numeric2mm5date(idate(nh))) if (idate(nh) == wrf_date(i)) then nh2wrf_file(nh) = wrf_file(i) i = i + 1 else nh2wrf_file(nh) = wrf_file(i-1) end if end do deallocate(wrf_file, wrf_date) ! Open 1st WRF file (NetCDF format) and read metadata if (chimverb > 2) print *, ' Opening 1st WRF file: ', trim(meteo_dir)//'/'//nh2wrf_file(0) call nc_get_struct(trim(meteo_dir)//'/'//nh2wrf_file(0), nc_WRF) call nc_get_meta(nc_WRF) do i=1, nc_WRF%ndim select case(trim(nc_WRF%dname(i))) case('west_east') nx_wrf = nc_WRF%dlen(i) case('south_north') ny_wrf = nc_WRF%dlen(i) case('bottom_top') nlev_wrf = nc_WRF%dlen(i) case('bottom_top_stag') nlev_F = nc_WRF%dlen(i) end select enddo numFrTimes = nc_WRF%nt nxx_DOT=nx_wrf+1 nyy_DOT=ny_wrf+1 ! Geometric calculations for vector interpolations if(.not.allocated(wrf_centers_3d))allocate(wrf_centers_3d(nx_wrf,ny_wrf,3)) if(.not.allocated(wrf_ulambda))allocate(wrf_ulambda(nx_wrf,ny_wrf,3)) if(.not.allocated(wrf_uphi))allocate(wrf_uphi(nx_wrf,ny_wrf,3)) if(.not.allocated(wrf_cosalpha))allocate(wrf_cosalpha(nx_wrf,ny_wrf)) if(.not.allocated(wrf_sinalpha))allocate(wrf_sinalpha(nx_wrf,ny_wrf)) wrf_cosalpha=nc_WRF%cosalpha wrf_sinalpha=nc_WRF%sinalpha call calculate_centers(nc_WRF%lon, nc_WRF%lat,wrf_centers_3d,wrf_ulambda,wrf_uphi,nx_wrf,ny_wrf) ! Are we using the same grid as WRF? iwrfgrid=test_wrf_grid(nc_wrf,nx_wrf,ny_wrf) if (iwrfgrid) then print*,'Runnig with offline meteorology on the same grid' else print*,'Runnig with offline meteorology on different grids' endif !-------------------------------- ! ONLINE else iwrfgrid = .true. ! Open WRFINPUT file (NetCDF format) and read metadata call nc_get_struct(fn_wrfinput, nc_WRFINPUT) call nc_get_meta(nc_WRFINPUT) do i=1, nc_WRFINPUT%ndim select case(trim(nc_WRFINPUT%dname(i))) case('bottom_top') nlev_wrf = nc_WRFINPUT%dlen(i) case('bottom_top_stag') nlev_F = nc_WRFINPUT%dlen(i) end select enddo nx_wrf = nzonal_entire_domain ny_wrf = nmerid_entire_domain nxx_DOT=nx_wrf+1 nyy_DOT=ny_wrf+1 ! Geometric calculations for vector interpolations if(.not.allocated(wrf_centers_3d))allocate(wrf_centers_3d(nx_wrf,ny_wrf,3)) if(.not.allocated(wrf_ulambda))allocate(wrf_ulambda(nx_wrf,ny_wrf,3)) if(.not.allocated(wrf_uphi))allocate(wrf_uphi(nx_wrf,ny_wrf,3)) if(.not.allocated(wrf_cosalpha))allocate(wrf_cosalpha(nx_wrf,ny_wrf)) if(.not.allocated(wrf_sinalpha))allocate(wrf_sinalpha(nx_wrf,ny_wrf)) wrf_cosalpha=nc_WRFINPUT%cosalpha wrf_sinalpha=nc_WRFINPUT%sinalpha wrf_centers_3d=chim_centers_3d wrf_ulambda=chim_ulambda wrf_uphi=chim_uphi ! Get ptop, sigmh amd sigmf in wrfinput file ! Caution : this assume that ptop is not a time varying parameter !!! allocate(sigmh_wrf(nlev_wrf)) allocate(sigmf_wrf(nlev_wrf+1)) ptop_wrf=nc_WRFINPUT%ptop sigmh_wrf=nc_WRFINPUT%znu sigmf_wrf=nc_WRFINPUT%znw endif !------------------------------ ! OFFLINE if (online == 0) then ! Get temporal increment in WRF files wrfDatestr = numeric2mm5date(nc_WRF%time(1)) read(wrfDatestr,'(I4,1x,I2,1x,I2,1x,I2,1x,I2,1x,I2)')& startYear,startMnth,startDay,startHour,startMin,startSec if(chimverb.ge.2)print*,'WRF file first date : '//wrfDatestr ! 2nd time slice to get the time step ! In the offline mode, wrf_dt is the timestep of the WRF input file, in hours wrfDatestr = numeric2mm5date(nc_WRF%time(2)) read(wrfDatestr,'(I4,1x,I2,1x,I2,1x,I2,1x,I2,1x,I2)')& startYear,startMnth,startDay,startHour1,startMin,startSec if (startHour1.ne.00) then wrf_dt=(startHour1-startHour)*60 else ihourcount=0 do it=1,numFrTimes wrfDatestr = numeric2mm5date(nc_WRF%time(it)) read(wrfDatestr,'(I4,1x,I2,1x,I2,1x,I2,1x,I2,1x,I2)')& & startYear,startMnth,startDay1,startHour1,startMin,startSec if (startDay1.eq.startDay.and.startHour1.eq.00) then ihourcount=ihourcount+1 endif enddo wrf_dt=60/ihourcount endif if(chimverb.ge.2)print*,'Meteo temporal increment: ',wrf_dt, ' minutes' ! number of time steps if(chimverb.ge.2)then print*,'Time steps in WRF file ',numFrTimes print*,'Time steps in CHIMERE ', nhours endif ! check end date/time of this file wrfDatestr = numeric2mm5date(nc_WRF%time(numFrTimes)) read(wrfDatestr,'(I4,1x,I2,1x,I2,1x,I2,1x,I2,1x,I2)')& endYear,endMnth,endDay,endHour,endMin,endSec ! Save Longitude/latitude of cells center call opfi(ifn_wrf,fn_mdom,'f','s') do ime=1,ny_wrf do izo=1,nx_wrf write(ifn_wrf,*) nc_WRF%lon(izo,ime), nc_WRF%lat(izo,ime) enddo enddo close(ifn_wrf) ! Save longitude/latitude of cell centers for U-wind call opfi(ifn_wrf,fn_mdom_u,'f','s') do ime=1,ny_wrf do izo=1,nxx_DOT write(ifn_wrf,*) nc_WRF%ulon(izo,ime), nc_WRF%ulat(izo,ime) enddo enddo close(ifn_wrf) ! Longitude/latitude of cell centers for V-wind call opfi(ifn_wrf,fn_mdom_v,'f','s') do ime=1,nyy_DOT do izo=1,nx_wrf write(ifn_wrf,*) nc_WRF%vlon(izo,ime), nc_WRF%vlat(izo,ime) enddo enddo close(ifn_wrf) !------------------------------ ! ONLINE else ! Save Longitude/latitude of cells center call opfi(ifn_wrf,fn_mdom,'f','s') do ime=1,ny_wrf do izo=1,nx_wrf write(ifn_wrf,*) xlong_entire_domain(izo,ime), xlati_entire_domain(izo,ime) enddo enddo close(ifn_wrf) endif ! Read weights and indexes for interpolation npts=nzonal_entire_domain*nmerid_entire_domain allocate(ixc(npts)) allocate(iyc(npts)) allocate(wxc(npts)) allocate(wyc(npts)) ! Cell centers for Mass fields if (.not.iwrfgrid) then if(chimverb.ge.2)print *, 'Calculation of weights for meteo interpolation ...' call prep_bilin(nx_wrf,ny_wrf,ixc,iyc,wxc,wyc,fn_mdom,fnigeo,npts) else i=0 do ime=1,nmerid_entire_domain do izo=1,nzonal_entire_domain i=i+1 ixc(i) = izo iyc(i) = ime enddo enddo ! These weights will not be used if iwrfgrid. If they were to be used the right values are 0. wxc=0. wyc=0. endif if (online == 0) then call nc_close(nc_WRF) else call nc_close(nc_WRFINPUT) endif end subroutine ini_prepwrf !------------------------------------------------------------------------------------ subroutine wrf_open_init_vars(fname, nc) ! Read WRF data in off-line mode only. ! Called by iniworker and read_wrf_file use chimere_common use pncvar implicit none character(len=*), intent(in) :: fname type(ncfile), intent(inout) :: nc integer :: ncstat, tmp logical :: is_present if (chimverb > 2 .and. rank == 0) print *,' Opening WRF file '//trim(fname)//'...' call pncvar_open(nc, trim(fname), NF90_NOWRITE, wrk_comm) call pncvar_get_times(nc) nx_wrf=findlen(nc,'west_east') ny_wrf=findlen(nc,'south_north') nlev_wrf=findlen(nc,'bottom_top') nvmet = findlen(nc,'bottom_top') call pncvar_get_var_info(nc, 'ZNU' , var_wrf_sigmh) call pncvar_get_var_info(nc, 'ZNW' , var_wrf_sigmf) call pncvar_get_var_info(nc, 'PSFC' , var_wrf_psfc) call pncvar_get_var_info(nc, 'SWDOWN', var_wrf_swrd) call pncvar_get_var_info(nc, 'GLW' , var_wrf_lwrd) call pncvar_get_var_info(nc, 'Q2' , var_wrf_r2m) call pncvar_get_var_info(nc, 'RAINC' , var_wrf_copc) call pncvar_get_var_info(nc, 'RAINNC', var_wrf_lspc) call pncvar_get_var_info(nc, 'LH' , var_wrf_slhf) call pncvar_get_var_info(nc, 'T2' , var_wrf_tem2) call pncvar_get_var_info(nc, 'U10' , var_wrf_u10) call pncvar_get_var_info(nc, 'V10' , var_wrf_v10) call pncvar_get_var_info(nc, 'SMOIS' , var_wrf_soim) call pncvar_get_var_info(nc, 'PB' , var_wrf_pref) call pncvar_get_var_info(nc, 'P' , var_wrf_prep) call pncvar_get_var_info(nc, 'U' , var_wrf_winx) call pncvar_get_var_info(nc, 'V' , var_wrf_winy) call pncvar_get_var_info(nc, 'W' , var_wrf_winw) call pncvar_get_var_info(nc, 'T' , var_wrf_thet) call pncvar_get_var_info(nc, 'QCLOUD', var_wrf_cliq) call pncvar_get_var_info(nc, 'QVAPOR', var_wrf_sphu) call pncvar_get_var_info(nc, 'PH' , var_wrf_ph) call pncvar_get_var_info(nc, 'PHB' , var_wrf_phb) ! Ice content not always contained in meteo files. Check it. cice_ok=1 call pncvar_inq_varid(nc%id,'QICE', tmp, var_is_present=is_present) if ( .not. is_present ) then cice_ok=0 else call pncvar_get_var_info(nc, 'QICE', var_wrf_cice) endif ! PHYD should be present in WRF files, we check it however. phyd_ok=1 call pncvar_inq_varid(nc%id,'P_HYD', tmp, var_is_present=is_present) if ( .not. is_present ) then phyd_ok=0 else call pncvar_get_var_info(nc, 'P_HYD', var_wrf_phyd) endif if(chimverb.ge.5.and.phyd_ok.eq.1.and.rank==0) then print*,'Hydrostatic pressure found. Use it for interpolations' else print*,'Warning : No hydrostatic pressure P_HYD in WRF outputs' endif ! Rain and snow should be present in WRF files, we check it however. rain_ok=1 call pncvar_inq_varid(nc%id,'QRAIN', tmp, var_is_present=is_present) if ( .not. is_present ) then rain_ok=0 else call pncvar_get_var_info(nc, 'QRAIN', var_wrf_rain) endif snow_ok=1 call pncvar_inq_varid(nc%id,'QSNOW', tmp, var_is_present=is_present) if ( .not. is_present ) then snow_ok=0 else call pncvar_get_var_info(nc, 'QSNOW', var_wrf_snow) endif fcld_ok=1 call pncvar_inq_varid(nc%id,'CLDFRA', tmp, var_is_present=is_present) if ( .not. is_present ) then fcld_ok=0 else call pncvar_get_var_info(nc, 'CLDFRA', var_wrf_fcld) endif ! Snow height not always contained in meteo files. Check it. snowh_ok=1 call pncvar_inq_varid(nc%id,'SNOWH', tmp, var_is_present=is_present) if ( .not. is_present ) then snowh_ok=0 else call pncvar_get_var_info(nc, 'SNOWH', var_wrf_snowh) endif ! Sea-ice content not always contained in meteo files. Check it. seaice_ok=1 call pncvar_inq_varid(nc%id,'SEAICE', tmp, var_is_present=is_present) if ( .not. is_present ) then seaice_ok=0 else call pncvar_get_var_info(nc, 'SEAICE', var_wrf_seaice) endif ! ustar not always contained in meteo files. Check it. usta_ok=1 call pncvar_inq_varid(nc%id,'UST', tmp, var_is_present=is_present) if ( .not. is_present ) then usta_ok=0 else call pncvar_get_var_info(nc, 'UST' , var_wrf_usta) endif ! sensible and latent heat fluxes flux_ok=1 call pncvar_inq_varid(nc%id,'HFX', tmp, var_is_present=is_present) if ( .not. is_present ) then flux_ok=0 else call pncvar_get_var_info(nc, 'HFX' , var_wrf_sshf) endif ! boundary layer height pblh_ok=1 call pncvar_inq_varid(nc%id,'PBLH', tmp, var_is_present=is_present) if ( .not. is_present ) then pblh_ok=0 else call pncvar_get_var_info(nc, 'PBLH' , var_wrf_hght) endif end subroutine wrf_open_init_vars !========================================================================================= subroutine read_wrf_file(itimeslot) ! Read WRF data in off-line mode only. use pncvar use chimere_common implicit none #define NCERR(lnum) if(ncstat/=NF90_NOERR) call nc_err(ncstat,lnum,__FILE__) integer,intent(in) :: itimeslot integer :: ncstat, i integer :: slot_wrf_file integer :: ptop_varid real(kind=iprec),dimension(:),allocatable :: tmp_vec integer,dimension(3) :: vec_count_3, vec_start_3 integer,dimension(4) :: vec_count_4, vec_start_4 logical :: new_file new_file = .true. if (itimeslot > 0) then if (trim(nh2wrf_file(itimeslot)) == trim(nh2wrf_file(itimeslot-1))) new_file = .false. endif if (new_file) then slot_wrf_file_prec = -1 ! Close previous file if (itimeslot > 0) call pncvar_close(wrf_nc) ! Open new file call wrf_open_init_vars(trim(meteo_dir)//'/'//trim(nh2wrf_file(itimeslot)), wrf_nc) ! Find 1st slot slot_wrf_file = -1 do i = 1, wrf_nc%nt if (wrf_nc%time(i) == idate(itimeslot)) then slot_wrf_file = i exit end if enddo if (slot_wrf_file == -1) then if (rank == 0) then print *, 'meteo-prepwrf : wrf file does not contain the requested time: ', idate(itimeslot) call exit1('Problematic file is '//trim(wrf_nc%name)) end if end if ini_slot_wrf_file = slot_wrf_file itimeslot0 = itimeslot else ! Pick next slot slot_wrf_file = 1 if (wrf_dt >= 60) then slot_wrf_file = ini_slot_wrf_file + (itimeslot - itimeslot0)/(wrf_dt/60) if (wrf_dt > 60 .and. mod(itimeslot*60, wrf_dt) > 0) slot_wrf_file = slot_wrf_file + 1 endif if (wrf_dt < 60) then slot_wrf_file = slot_wrf_file_prec + (60/wrf_dt) endif end if if (rank == 0 .and. chimverb > 5) print *, ' -- itimeslot =', itimeslot, & & ' slot_wrf_file =', slot_wrf_file, ' ini_slot_wrf_file =', ini_slot_wrf_file ! If the current slot in wrf file has not changed it is not necesary to re-read it if (itimeslot /= 0) then if(slot_wrf_file_prec == slot_wrf_file)return ! When itimeslot=0 slot_wrf_file_prec is not initialized endif slot_wrf_file_prec = slot_wrf_file ! Read WRF top pressure ! Caution : this may be a time varying parameter !!! call pncvar_inq_varid(wrf_nc%id, 'P_TOP', ptop_varid) allocate(tmp_vec(1)) call pncvar_get_vara_double(wrf_nc%id, ptop_varid, (/1/), (/1/), tmp_vec) ptop_wrf=tmp_vec(1) deallocate(tmp_vec) ! Get variables call pncvar_get_var(domains,var_wrf_sigmh, sigmh_wrf(1:nlev_wrf), & & rank,(/nlev_wrf,1/),(/1,slot_wrf_file/)) call pncvar_get_var(domains,var_wrf_sigmf, sigmf_wrf(1:nlev_wrf+1), & & rank,(/nlev_wrf+1,1/),(/1,slot_wrf_file/)) vec_count_3 = (/nzcount_wrf,nmcount_wrf,1/) vec_start_3 = (/izstart_chim,imstart_chim,slot_wrf_file/) call pncvar_get_var(domains,var_wrf_psfc, & & psfc_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) call pncvar_get_var(domains,var_wrf_swrd, & & swrd_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) call pncvar_get_var(domains,var_wrf_lwrd, & & lwrd_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) call pncvar_get_var(domains,var_wrf_r2m, & & r2m_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) call pncvar_get_var(domains,var_wrf_copc, & & copc_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) call pncvar_get_var(domains,var_wrf_lspc, & & lspc_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) if (flux_ok==1) call pncvar_get_var(domains,var_wrf_sshf, & & sshf_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) call pncvar_get_var(domains,var_wrf_slhf, & & slhf_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) if (usta_ok==1) call pncvar_get_var(domains,var_wrf_usta, & & usta_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) call pncvar_get_var(domains,var_wrf_tem2, & & tem2_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) if (pblh_ok==1) call pncvar_get_var(domains,var_wrf_hght, & & hght_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) if(snowh_ok==1)call pncvar_get_var(domains,var_wrf_snowh, & & snowh_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) if(seaice_ok==1)call pncvar_get_var(domains,var_wrf_seaice, & & seaice_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) call pncvar_get_var(domains,var_wrf_u10, & & u10_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) call pncvar_get_var(domains,var_wrf_v10, & & v10_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_3,vec_start_3,.false.) vec_count_4 = (/nzcount_wrf,nmcount_wrf,nlev_wrf,1/) vec_start_4 = (/izstart_chim,imstart_chim,1,slot_wrf_file/) call pncvar_get_var(domains,var_wrf_pref, & & pref_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,1:nlev_wrf), & & rank,vec_count_4,vec_start_4,.false.) call pncvar_get_var(domains,var_wrf_prep, & & prep_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,1:nlev_wrf), & & rank,vec_count_4,vec_start_4,.false.) if (phyd_ok.eq.1)then call pncvar_get_var(domains,var_wrf_phyd, & & phyd_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,1:nlev_wrf), & & rank,vec_count_4,vec_start_4,.false.) endif call pncvar_get_var(domains,var_wrf_thet, & & thet_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,1:nlev_wrf), & & rank,vec_count_4,vec_start_4,.false.) call pncvar_get_var(domains,var_wrf_cliq, & & cliq_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,1:nlev_wrf), & & rank,vec_count_4,vec_start_4,.false.) call pncvar_get_var(domains,var_wrf_rain, & & rain_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,1:nlev_wrf), & & rank,vec_count_4,vec_start_4,.false.) if(snow_ok==1) then call pncvar_get_var(domains,var_wrf_snow, & & snow_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,1:nlev_wrf), & & rank,vec_count_4,vec_start_4,.false.) else snow_wrf=dzero endif if(fcld_ok.eq.1)then call pncvar_get_var(domains,var_wrf_fcld, & & fcld_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,1:nlev_wrf), & & rank,vec_count_4,vec_start_4,.false.) else fcld_wrf=dzero endif call pncvar_get_var(domains,var_wrf_sphu, & & sphu_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,1:nlev_wrf), & & rank,vec_count_4,vec_start_4,.false.) if(cice_ok==1) call pncvar_get_var(domains,var_wrf_cice, & & cice_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,1:nlev_wrf), & & rank,vec_count_4,vec_start_4,.false.) vec_count_4 = (/nzcount_wrf,nmcount_wrf,nlev_wrf+1,1/) vec_start_4 = (/izstart_chim,imstart_chim,1,slot_wrf_file/) call pncvar_get_var(domains,var_wrf_ph, & & ph_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,1:nlev_wrf+1), & & rank,vec_count_4,vec_start_4,.false.) call pncvar_get_var(domains,var_wrf_phb, & & phb_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,1:nlev_wrf+1), & & rank,vec_count_4,vec_start_4,.false.) vec_count_4 = (/nzcount_wrf,nmcount_wrf,1,1/) vec_start_4 = (/izstart_chim,imstart_chim,1,slot_wrf_file/) call pncvar_get_var(domains,var_wrf_soim, & & soim_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim), & & rank,vec_count_4,vec_start_4,.false.) call pncvar_get_var(domains,var_wrf_winx, & & winx_wrf(izstart_chim:izend_chim+1,imstart_chim:imend_chim,1:nlev_wrf), & & rank,(/nzcount_wrf+1,nmcount_wrf,nlev_wrf,1/), & & (/izstart_chim,imstart_chim,1,slot_wrf_file/),.false.) call pncvar_get_var(domains,var_wrf_winy, & & winy_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim+1,1:nlev_wrf), & & rank,(/nzcount_wrf,nmcount_wrf+1,nlev_wrf,1/), & & (/izstart_chim,imstart_chim,1,slot_wrf_file/),.false.) call pncvar_get_var(domains,var_wrf_winw, & & winw_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,1:nlev_wrf), & & rank,(/nzcount_wrf,nmcount_wrf,nlev_wrf,1/), & & (/izstart_chim,imstart_chim,1,slot_wrf_file/),.false.) end subroutine read_wrf_file !================================================================================================ subroutine convert_wrf_fields(ks,itimeslot) ! This routine is used both in off-line and on-line mode use chimere_common use interp_tools implicit none #define NCERR(lnum) if(ncstat/=NF90_NOERR) call nc_err(ncstat,lnum,__FILE__) integer,intent(in) :: itimeslot,ks integer :: i,j,k, k2, wrf_dt_h real(kind=iprec) :: z1,z2,orog,dz,logfull,loghalf real(kind=iprec),dimension(:,:),allocatable :: pstar real(kind=iprec),dimension(:,:,:),allocatable :: alti_wrf,temp_wrf,winz_wrf,winm_wrf real(kind=iprec),dimension(:,:),allocatable :: u10mz_wrf,v10mm_wrf ! Earth-coordinate (zonal and meridional) 10m wind components ! T00 here is NOT the same as that found in WRF output TOO variable !! ! See the wrf-users mail list answers of 12/06/2013 and WRF-ARW documentation ! => Use always the constant +300 to get the total potential temperature t00 = 300. ! cice_wrf, cliq_wrf et sphu_wrf shall be positive where (cliq_wrf .lt. 0) cliq_wrf = 0. where (snowh_wrf .lt. 0) snowh_wrf = 0. where (rain_wrf .lt. 0) rain_wrf = 0. where (snow_wrf .lt. 0) snow_wrf = 0. where (fcld_wrf .lt. 0) fcld_wrf = 0. if(cice_ok==1) where (cice_wrf .lt. 0) cice_wrf = 0. ! if seaice is present, impose 10 cm of snow on top of ocean to have the correct albedo if((seaice_ok==1).and.(snowh_ok==1)) where (seaice_wrf .gt. 0.5) snowh_wrf = 0.1 ! too small specific humidities are not allowed due to a log(sphu) in Chimere. where (sphu_wrf .lt. min_sphu_wrf) sphu_wrf = min_sphu_wrf ! All WRF variables, at the exception of Vertical Velocity and Geopotential ! are defined at half-sigma levels. ! Let us compute the altitudes of these half-sigma levels. ! altitude of middle of layers (output) if (.not. allocated(alti_wrf)) allocate(alti_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,nlev_wrf)) if (.not. allocated(pstar)) allocate(pstar(izstart_chim:izend_chim,imstart_chim:imend_chim)) if ( ptop_wrf.lt.0.1.or.ptop_wrf.ne.ptop_wrf) then print*, __FILE__//': ERROR: PTOP = ',ptop_wrf, 'Bye !' call exit1('Exiting') endif do j=imstart_chim,imend_chim do i=izstart_chim,izend_chim pstar(i,j) = psfc_wrf(i,j) - ptop_wrf ! orography height orog = (ph_wrf(i,j,1) + phb_wrf(i,j,1)) /gravit do k=1,nlev_wrf logfull=log((pstar(i,j)*sigmf_wrf(k)+ptop_wrf)/(pstar(i,j)*sigmf_wrf(k+1)+ptop_wrf)) loghalf=log((pstar(i,j)*sigmf_wrf(k)+ptop_wrf)/(pstar(i,j)*sigmh_wrf(k) +ptop_wrf)) z2 = (ph_wrf(i,j,k+1) + phb_wrf(i,j,k+1)) /gravit z1 = (ph_wrf(i,j,k) + phb_wrf(i,j,k)) /gravit if(logfull.eq.0)then print*,'Warning: logfull is zero in rank',rank,' and cell',i,j,k dz = (z2 - z1) else dz = (z2 - z1) * loghalf / logfull endif ! altitude of half-sigma levels alti_wrf(i,j,k) = z1 + dz - orog if(alti_wrf(i,j,k).lt.0)then print*,'WARNING: Negative WRF altitude in',i,j,k,alti_wrf(i,j,k) endif if(alti_wrf(i,j,k) /= alti_wrf(i,j,k))then print*,'-----------------------------------------' print*,'alti_wrf problem in meteo-prepwrf rank',rank print*,'NaN values received/read from the met model' print*,'If this an online simulation, try reducing WRF timestep' print*,'Some debug information:' print*,'i j k',i,j,k print*,'ph_wrf k',ph_wrf(i,j,k),' phb_wrf k',phb_wrf(i,j,k) print*,'ph_wrf k+1',ph_wrf(i,j,k+1),' phb_wrf k+1',phb_wrf(i,j,k+1) print*,'z1=',z1,' z2=',z2,' dz=',dz print*,'-----------------------------------------' call exit1('NaN values for altitude. STOPPING') endif enddo enddo enddo if(allocated(pstar))deallocate(pstar) ! Pressure and temperature if (.not.allocated(temp_wrf))allocate(temp_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,nlev_wrf)) ! lm: Check values of P and T: to remove later do j=imstart_chim,imend_chim do i=izstart_chim,izend_chim do k=1,nlev_wrf pres_wrf(i,j,k) = pref_wrf(i,j,k) + prep_wrf(i,j,k) temp_wrf(i,j,k) = (thet_wrf(i,j,k)+t00) * (pres_wrf(i,j,k)/p00)**xkappa if(temp_wrf(i,j,k).le.100.or.temp_wrf(i,j,k).ge.350)then print*,'meteo-prepwrf warning: extreme temperature',i,j,k print*,' temp=theta+t00 ',temp_wrf(i,j,k),thet_wrf(i,j,k),t00 print*,' pres=pref+prep ',pres_wrf(i,j,k),pref_wrf(i,j,k),prep_wrf(i,j,k) print*,' vertical profile:' do k2=1,nlev_wrf print*,k2,prep_wrf(i,j,k2),pres_wrf(i,j,k2),temp_wrf(i,j,k2),thet_wrf(i,j,k2) enddo endif enddo enddo enddo ! All pertinent variables have been stored into buffers. ! It is time to perform some calculations, and to write the result to the output file. ! for all useful time itimeslots ! Do U and V rotation for winds allocate(winz_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,nlev_wrf)) allocate(winm_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim,nlev_wrf)) call do_uv_rot_sm(winx_wrf,winy_wrf,winz_wrf,winm_wrf, & & wrf_cosalpha(izstart_chim:izend_chim,imstart_chim:imend_chim), & & wrf_sinalpha(izstart_chim:izend_chim,imstart_chim:imend_chim)) ! At this stage of the processing, we have the variables for a time itimeslot ! stored into memory, in the WRF grid coordinates system. ! Now we interpolate each variable onto the chimere grid coordinates system. ! Horizontal interpolation call wrf2chim_2d(psfc_wrf,m_psfc(1:nzonal,1:nmerid,:),ks) call wrf2chim_2d(swrd_wrf,m_swrd(1:nzonal,1:nmerid,:),ks) call wrf2chim_2d(lwrd_wrf,m_lwrd(1:nzonal,1:nmerid,:),ks) call wrf2chim_2d(tem2_wrf,m_tem2(1:nzonal,1:nmerid,:),ks) call wrf2chim_2d(soim_wrf,m_soim(1:nzonal,1:nmerid,:),ks) ! Precips must be calculated for each time slot_wrf_file since the beginning of the WRF file. ! WRF provides accumulated fields, and chimere requires differences. ! precip at t=0 always=0 will not be used in chimere anyway: interp of precip ! is done with cumulated precip up to the next hour in locvalues if (.not. allocated(wrf_cop0)) allocate(wrf_cop0(izstart_chim:izend_chim,imstart_chim:imend_chim)) if (.not. allocated(wrf_lsp0)) allocate(wrf_lsp0(izstart_chim:izend_chim,imstart_chim:imend_chim)) ! For cumulated variables, need to initialize with the first field, not zero if (itimeslot == 0) then wrf_lsp0=lspc_wrf wrf_cop0=copc_wrf endif call wrf2chim_2d(lspc_wrf - wrf_lsp0,m_lspc(1:nzonal,1:nmerid,:),ks) call wrf2chim_2d(copc_wrf - wrf_cop0,m_copc(1:nzonal,1:nmerid,:),ks) ! 'if' test has to be splitted due to the non-sense of wrf_dt in online mode. if(ks.eq.2)then if(online.eq.1)then wrf_lsp0(izstart_chim:izend_chim,imstart_chim:imend_chim) = & & lspc_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim) wrf_cop0(izstart_chim:izend_chim,imstart_chim:imend_chim) = & & copc_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim) else if(mod(itimeslot*60,wrf_dt).eq.0)then wrf_lsp0(izstart_chim:izend_chim,imstart_chim:imend_chim) = & & lspc_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim) wrf_cop0(izstart_chim:izend_chim,imstart_chim:imend_chim) = & & copc_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim) endif endif endif ! u10m and v10m ! rotation to Earth-system coordinates for u10m and v10m allocate(u10mz_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim)) allocate(v10mm_wrf(izstart_chim:izend_chim,imstart_chim:imend_chim)) call do_uv_rot_2d(u10_wrf,v10_wrf,u10mz_wrf,v10mm_wrf, & & wrf_cosalpha(izstart_chim:izend_chim,imstart_chim:imend_chim), & & wrf_sinalpha(izstart_chim:izend_chim,imstart_chim:imend_chim)) ! Interpolate 10m wind vectors onto the CHIMERE grid call wrf2chim_vec_2d(u10mz_wrf,v10mm_wrf,m_u10m(1:nzonal,1:nmerid,:),m_v10m(1:nzonal,1:nmerid,:),ks) deallocate(u10mz_wrf) deallocate(v10mm_wrf) ! snowh if(snowh_ok==1)then call wrf2chim_2d(snowh_wrf,m_snowh(1:nzonal,1:nmerid,:),ks) else m_snowh(1:nzonal,1:nmerid,:)=0. endif ! seaice if(seaice_ok==1)then call wrf2chim_2d(seaice_wrf,m_seaice(1:nzonal,1:nmerid,:),ks) else m_seaice(1:nzonal,1:nmerid,:)=0. endif ! sshf if (flux_ok==1) then call wrf2chim_2d(sshf_wrf,m_sshf(1:nzonal,1:nmerid,:),ks) else m_sshf(1:nzonal,1:nmerid,:) = 0. endif ! slhf call wrf2chim_2d(slhf_wrf,m_slhf(1:nzonal,1:nmerid,:),ks) ! usta if (usta_ok==1) then call wrf2chim_2d(usta_wrf,m_usta(1:nzonal,1:nmerid,:),ks) else m_usta(1:nzonal,1:nmerid,ks) = 0. endif ! hght if (pblh_ok == 1) then call wrf2chim_2d(hght_wrf,m_hght(1:nzonal,1:nmerid,:),ks) else m_hght(1:nzonal,1:nmerid,ks) = 0. endif ! 3d variables call wrf2chim_3d(temp_wrf,m_temp(1:nzonal,1:nmerid,:,:),ks) call wrf2chim_3d(cliq_wrf,m_cliq(1:nzonal,1:nmerid,:,:),ks) call wrf2chim_3d(rain_wrf,m_rain(1:nzonal,1:nmerid,:,:),ks) call wrf2chim_3d(fcld_wrf,m_fcld(1:nzonal,1:nmerid,:,:),ks) call wrf2chim_3d(sphu_wrf,m_sphu(1:nzonal,1:nmerid,:,:),ks) call wrf2chim_3d(pres_wrf,m_pres(1:nzonal,1:nmerid,:,:),ks) call wrf2chim_3d(alti_wrf,m_alti(1:nzonal,1:nmerid,:,:),ks) if(phyd_ok.eq.1)then call wrf2chim_3d(phyd_wrf,m_phyd(1:nzonal,1:nmerid,:,:),ks) endif if (.not.(is_diagwinw)) call wrf2chim_3d(winw_wrf,m_winw(1:nzonal,1:nmerid,:,:),ks) ! m_winz, m_winm call wrf2chim_vec_3d(winz_wrf,winm_wrf,m_winz(1:nzonal,1:nmerid,:,:),m_winm(1:nzonal,1:nmerid,:,:),ks) if(cice_ok==1) then call wrf2chim_3d(cice_wrf,m_cice(1:nzonal,1:nmerid,:,:),ks) else m_cice(1:nzonal,1:nmerid,:,ks) = 0. endif if(snow_ok==1) then call wrf2chim_3d(snow_wrf,m_snow(1:nzonal,1:nmerid,:,:),ks) else m_snow(1:nzonal,1:nmerid,:,ks) = dzero endif ! Save previous time step for temporal interpolation if (itimeslot == 0) then m_psfc_prec(1:nzonal,1:nmerid) = m_psfc(1:nzonal,1:nmerid,ks) m_swrd_prec(1:nzonal,1:nmerid) = m_swrd(1:nzonal,1:nmerid,ks) m_lwrd_prec(1:nzonal,1:nmerid) = m_lwrd(1:nzonal,1:nmerid,ks) m_tem2_prec(1:nzonal,1:nmerid) = m_tem2(1:nzonal,1:nmerid,ks) m_soim_prec(1:nzonal,1:nmerid) = m_soim(1:nzonal,1:nmerid,ks) m_u10m_prec(1:nzonal,1:nmerid) = m_u10m(1:nzonal,1:nmerid,ks) m_v10m_prec(1:nzonal,1:nmerid) = m_v10m(1:nzonal,1:nmerid,ks) if(snowh_ok==1)m_snowh_prec(1:nzonal,1:nmerid) = m_snowh(1:nzonal,1:nmerid,ks) m_seaice_prec(1:nzonal,1:nmerid) = m_seaice(1:nzonal,1:nmerid,ks) m_sshf_prec(1:nzonal,1:nmerid) = m_sshf(1:nzonal,1:nmerid,ks) m_slhf_prec(1:nzonal,1:nmerid) = m_slhf(1:nzonal,1:nmerid,ks) m_usta_prec(1:nzonal,1:nmerid) = m_usta(1:nzonal,1:nmerid,ks) m_hght_prec(1:nzonal,1:nmerid) = m_hght(1:nzonal,1:nmerid,ks) m_temp_prec(1:nzonal,1:nmerid,1:nlev_wrf) = m_temp(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_cliq_prec(1:nzonal,1:nmerid,1:nlev_wrf) = m_cliq(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_rain_prec(1:nzonal,1:nmerid,1:nlev_wrf) = m_rain(1:nzonal,1:nmerid,1:nlev_wrf,ks) if(snow_ok==1)m_snow_prec(1:nzonal,1:nmerid,1:nlev_wrf) = m_snow(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_fcld_prec(1:nzonal,1:nmerid,1:nlev_wrf) = m_fcld(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_sphu_prec(1:nzonal,1:nmerid,1:nlev_wrf) = m_sphu(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_pres_prec(1:nzonal,1:nmerid,1:nlev_wrf) = m_pres(1:nzonal,1:nmerid,1:nlev_wrf,ks) if(phyd_ok==1)m_phyd_prec(1:nzonal,1:nmerid,1:nlev_wrf) = m_phyd(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_alti_prec(1:nzonal,1:nmerid,1:nlev_wrf) = m_alti(1:nzonal,1:nmerid,1:nlev_wrf,ks) if(cice_ok==1)m_cice_prec(1:nzonal,1:nmerid,1:nlev_wrf) = m_cice(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_winz_prec(1:nzonal,1:nmerid,1:nlev_wrf) = m_winz(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_winm_prec(1:nzonal,1:nmerid,1:nlev_wrf) = m_winm(1:nzonal,1:nmerid,1:nlev_wrf,ks) endif ! Temporal interpolation, only in offline mode if(online.eq.0)then if (ks == 2 .and. wrf_dt > 60 ) then wrf_dt_h=wrf_dt/60 dt=mod(itimeslot-1,wrf_dt_h)+1 m_psfc(1:nzonal,1:nmerid,ks)=(m_psfc_prec(1:nzonal,1:nmerid)*(wrf_dt_h-dt) & & +m_psfc(1:nzonal,1:nmerid,ks)*dt)/wrf_dt_h m_swrd(1:nzonal,1:nmerid,ks)=(m_swrd_prec(1:nzonal,1:nmerid)*(wrf_dt_h-dt) & & +m_swrd(1:nzonal,1:nmerid,ks)*dt)/wrf_dt_h m_lwrd(1:nzonal,1:nmerid,ks)=(m_lwrd_prec(1:nzonal,1:nmerid)*(wrf_dt_h-dt) & & +m_lwrd(1:nzonal,1:nmerid,ks)*dt)/wrf_dt_h m_tem2(1:nzonal,1:nmerid,ks)=(m_tem2_prec(1:nzonal,1:nmerid)*(wrf_dt_h-dt) & & +m_tem2(1:nzonal,1:nmerid,ks)*dt)/wrf_dt_h m_soim(1:nzonal,1:nmerid,ks)=(m_soim_prec(1:nzonal,1:nmerid)*(wrf_dt_h-dt) & & +m_soim(1:nzonal,1:nmerid,ks)*dt)/wrf_dt_h m_u10m(1:nzonal,1:nmerid,ks)=(m_u10m_prec(1:nzonal,1:nmerid)*(wrf_dt_h-dt) & & +m_u10m(1:nzonal,1:nmerid,ks)*dt)/wrf_dt_h m_v10m(1:nzonal,1:nmerid,ks)=(m_v10m_prec(1:nzonal,1:nmerid)*(wrf_dt_h-dt) & & +m_v10m(1:nzonal,1:nmerid,ks)*dt)/wrf_dt_h if(snowh_ok==1)m_snowh(1:nzonal,1:nmerid,ks)=(m_snowh_prec(1:nzonal,1:nmerid)*(wrf_dt_h-dt) & & +m_snowh(1:nzonal,1:nmerid,ks)*dt)/wrf_dt_h if(seaice_ok==1)m_seaice(1:nzonal,1:nmerid,ks)=(m_seaice_prec(1:nzonal,1:nmerid)*(wrf_dt_h-dt) & & +m_seaice(1:nzonal,1:nmerid,ks)*dt)/wrf_dt_h m_sshf(1:nzonal,1:nmerid,ks)=(m_sshf_prec(1:nzonal,1:nmerid)*(wrf_dt_h-dt) & & +m_sshf(1:nzonal,1:nmerid,ks)*dt)/wrf_dt_h m_slhf(1:nzonal,1:nmerid,ks)=(m_slhf_prec(1:nzonal,1:nmerid)*(wrf_dt_h-dt) & & +m_slhf(1:nzonal,1:nmerid,ks)*dt)/wrf_dt_h m_usta(1:nzonal,1:nmerid,ks)=(m_usta_prec(1:nzonal,1:nmerid)*(wrf_dt_h-dt) & & +m_usta(1:nzonal,1:nmerid,ks)*dt)/wrf_dt_h m_hght(1:nzonal,1:nmerid,ks)=(m_hght_prec(1:nzonal,1:nmerid)*(wrf_dt_h-dt) & & +m_hght(1:nzonal,1:nmerid,ks)*dt)/wrf_dt_h m_temp(1:nzonal,1:nmerid,1:nlev_wrf,ks)=(m_temp_prec(1:nzonal,1:nmerid,1:nlev_wrf)*(wrf_dt_h-dt)& & +m_temp(1:nzonal,1:nmerid,1:nlev_wrf,ks)*dt)/wrf_dt_h m_cliq(1:nzonal,1:nmerid,1:nlev_wrf,ks)=(m_cliq_prec(1:nzonal,1:nmerid,1:nlev_wrf)*(wrf_dt_h-dt)& & +m_cliq(1:nzonal,1:nmerid,1:nlev_wrf,ks)*dt)/wrf_dt_h m_rain(1:nzonal,1:nmerid,1:nlev_wrf,ks)=(m_rain_prec(1:nzonal,1:nmerid,1:nlev_wrf)*(wrf_dt_h-dt)& & +m_rain(1:nzonal,1:nmerid,1:nlev_wrf,ks)*dt)/wrf_dt_h if(snow_ok==1)m_snow(1:nzonal,1:nmerid,1:nlev_wrf,ks)=(m_snow_prec(1:nzonal,1:nmerid,1:nlev_wrf)*(wrf_dt_h-dt)& & +m_snow(1:nzonal,1:nmerid,1:nlev_wrf,ks)*dt)/wrf_dt_h m_fcld(1:nzonal,1:nmerid,1:nlev_wrf,ks)=(m_fcld_prec(1:nzonal,1:nmerid,1:nlev_wrf)*(wrf_dt_h-dt)& & +m_fcld(1:nzonal,1:nmerid,1:nlev_wrf,ks)*dt)/wrf_dt_h m_sphu(1:nzonal,1:nmerid,1:nlev_wrf,ks)=(m_sphu_prec(1:nzonal,1:nmerid,1:nlev_wrf)*(wrf_dt_h-dt)& & +m_sphu(1:nzonal,1:nmerid,1:nlev_wrf,ks)*dt)/wrf_dt_h m_pres(1:nzonal,1:nmerid,1:nlev_wrf,ks)=(m_pres_prec(1:nzonal,1:nmerid,1:nlev_wrf)*(wrf_dt_h-dt)& & +m_pres(1:nzonal,1:nmerid,1:nlev_wrf,ks)*dt)/wrf_dt_h if(phyd_ok==1)m_phyd(1:nzonal,1:nmerid,1:nlev_wrf,ks)=(m_phyd_prec(1:nzonal,1:nmerid,1:nlev_wrf)*(wrf_dt_h-dt)& & +m_phyd(1:nzonal,1:nmerid,1:nlev_wrf,ks)*dt)/wrf_dt_h m_alti(1:nzonal,1:nmerid,1:nlev_wrf,ks)=(m_alti_prec(1:nzonal,1:nmerid,1:nlev_wrf)*(wrf_dt_h-dt)& & +m_alti(1:nzonal,1:nmerid,1:nlev_wrf,ks)*dt)/wrf_dt_h m_cice(1:nzonal,1:nmerid,1:nlev_wrf,ks)=(m_cice_prec(1:nzonal,1:nmerid,1:nlev_wrf)*(wrf_dt_h-dt)& & +m_cice(1:nzonal,1:nmerid,1:nlev_wrf,ks)*dt)/wrf_dt_h m_winz(1:nzonal,1:nmerid,1:nlev_wrf,ks)=(m_winz_prec(1:nzonal,1:nmerid,1:nlev_wrf)*(wrf_dt_h-dt)& & +m_winz(1:nzonal,1:nmerid,1:nlev_wrf,ks)*dt)/wrf_dt_h m_winm(1:nzonal,1:nmerid,1:nlev_wrf,ks)=(m_winm_prec(1:nzonal,1:nmerid,1:nlev_wrf)*(wrf_dt_h-dt)& & +m_winm(1:nzonal,1:nmerid,1:nlev_wrf,ks)*dt)/wrf_dt_h ! special case for precip: accumulated since last output in WRF ! if the output is not hourly: convert to mm/hr m_lspc(1:nzonal,1:nmerid,ks)=m_lspc(1:nzonal,1:nmerid,ks)/wrf_dt_h m_copc(1:nzonal,1:nmerid,ks)=m_copc(1:nzonal,1:nmerid,ks)/wrf_dt_h ! Previous current time step data only if a new slot will be read at the new time step. if (mod(itimeslot,wrf_dt_h) == 0) then m_psfc_prec(1:nzonal,1:nmerid) = m_psfc(1:nzonal,1:nmerid,ks) m_swrd_prec(1:nzonal,1:nmerid) = m_swrd(1:nzonal,1:nmerid,ks) m_lwrd_prec(1:nzonal,1:nmerid) = m_lwrd(1:nzonal,1:nmerid,ks) m_tem2_prec(1:nzonal,1:nmerid) = m_tem2(1:nzonal,1:nmerid,ks) m_soim_prec(1:nzonal,1:nmerid) = m_soim(1:nzonal,1:nmerid,ks) m_u10m_prec(1:nzonal,1:nmerid) = m_u10m(1:nzonal,1:nmerid,ks) m_v10m_prec(1:nzonal,1:nmerid) = m_v10m(1:nzonal,1:nmerid,ks) m_snowh_prec(1:nzonal,1:nmerid) = m_snowh(1:nzonal,1:nmerid,ks) if(seaice_ok==1)m_seaice_prec(1:nzonal,1:nmerid) = m_seaice(1:nzonal,1:nmerid,ks) m_sshf_prec(1:nzonal,1:nmerid) = m_sshf(1:nzonal,1:nmerid,ks) m_slhf_prec(1:nzonal,1:nmerid) = m_slhf(1:nzonal,1:nmerid,ks) m_usta_prec(1:nzonal,1:nmerid) = m_usta(1:nzonal,1:nmerid,ks) m_hght_prec(1:nzonal,1:nmerid) = m_hght(1:nzonal,1:nmerid,ks) m_temp_prec(1:nzonal,1:nmerid, 1:nlev_wrf) = m_temp(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_cliq_prec(1:nzonal,1:nmerid, 1:nlev_wrf) = m_cliq(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_rain_prec(1:nzonal,1:nmerid, 1:nlev_wrf) = m_rain(1:nzonal,1:nmerid,1:nlev_wrf,ks) if(snow_ok==1)m_snow_prec(1:nzonal,1:nmerid, 1:nlev_wrf) = m_snow(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_fcld_prec(1:nzonal,1:nmerid, 1:nlev_wrf) = m_fcld(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_sphu_prec(1:nzonal,1:nmerid, 1:nlev_wrf) = m_sphu(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_pres_prec(1:nzonal,1:nmerid, 1:nlev_wrf) = m_pres(1:nzonal,1:nmerid,1:nlev_wrf,ks) if(phyd_ok==1)m_phyd_prec(1:nzonal,1:nmerid, 1:nlev_wrf) = m_phyd(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_alti_prec(1:nzonal,1:nmerid, 1:nlev_wrf) = m_alti(1:nzonal,1:nmerid,1:nlev_wrf,ks) if(cice_ok==1)m_cice_prec(1:nzonal,1:nmerid, 1:nlev_wrf) = m_cice(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_winz_prec(1:nzonal,1:nmerid, 1:nlev_wrf) = m_winz(1:nzonal,1:nmerid,1:nlev_wrf,ks) m_winm_prec(1:nzonal,1:nmerid, 1:nlev_wrf) = m_winm(1:nzonal,1:nmerid,1:nlev_wrf,ks) endif endif endif end subroutine convert_wrf_fields !----------------------------------------------------------------------------------------------- subroutine wrf2chim_2d(wrf_field,chim_field,ks) use chimere_common implicit none integer,intent(in) :: ks real(kind=iprec),intent(in) :: wrf_field(izend_chim-izstart_chim+1,imend_chim-imstart_chim+1) real(kind=iprec),intent(inout) :: chim_field(nzonal,nmerid,2) integer :: index_zonal,index_merid,np,izo,ime do ime=1,nmerid do izo=1,nzonal np=domains%dom(rank+1)%izstart+izo+(domains%dom(rank+1)%imstart+ime-2)*nzonal_entire_domain-1 index_zonal=ixc(np)-izstart_chim+1 index_merid=iyc(np)-imstart_chim+1 if(rank==0)then if(index_zonal.gt.izend_chim-izstart_chim+1)print*,'index_zonal out of bound' if(index_merid.gt.imend_chim-imstart_chim+1)print*,'index_merid out of bound' endif if ((index_zonal.ne.izend_chim-izstart_chim+1).and.(index_merid.ne.imend_chim-imstart_chim+1).and. (.not.iwrfgrid)) then chim_field(izo,ime,ks) = (1.-wxc(np))*(1.-wyc(np))*wrf_field(index_zonal,index_merid) & & + wxc(np) *(1.-wyc(np))*wrf_field(index_zonal+1,index_merid) & & + (1.-wxc(np))* wyc(np) *wrf_field(index_zonal, index_merid+1) & & + wxc(np) * wyc(np) *wrf_field(index_zonal+1,index_merid+1) else chim_field(izo,ime,ks) = wrf_field(index_zonal,index_merid) endif enddo enddo end subroutine wrf2chim_2d !----------------------------------------------------------------------------------------------- subroutine wrf2chim_3d(wrf_field,chim_field,ks) !chim_field inout because we want to respect slot 1 when filling slot 2 use chimere_common implicit none integer,intent(in) :: ks real(kind=iprec),intent(in) :: wrf_field(izend_chim-izstart_chim+1,imend_chim-imstart_chim+1,nlev_wrf) real(kind=iprec),intent(inout) :: chim_field(nzonal,nmerid,nlev_wrf,2) integer :: iz do iz=1,nlev_wrf call wrf2chim_2d(wrf_field(:,:,iz),chim_field(:,:,iz,:),ks) enddo end subroutine wrf2chim_3d !----------------------------------------------------------------------------------------------- subroutine wrf2chim_vec_2d(wrf_u,wrf_v,chim_u,chim_v,ks) ! inout because we want to respect slot 1 when filling slot 2 use interp_tools, only : scalar use chimere_common implicit none integer,intent(in) :: ks real(kind=iprec),intent(in) :: wrf_u(izend_chim-izstart_chim+1,imend_chim-imstart_chim+1) real(kind=iprec),intent(in) :: wrf_v(izend_chim-izstart_chim+1,imend_chim-imstart_chim+1) real(kind=iprec),intent(inout) :: chim_u(nzonal,nmerid,2),chim_v(nzonal,nmerid,2) integer :: index_zonal,index_merid,np,ix,iy,izo,ime real(kind=iprec) :: w3d(izend_chim-izstart_chim+1,imend_chim-imstart_chim+1,3) real(kind=iprec) :: w3dout(nzonal,nmerid,3) do ix=1,izend_chim-izstart_chim+1 do iy=1, imend_chim-imstart_chim+1 w3d(ix,iy,:)= wrf_u(ix,iy) *wrf_ulambda(ix+izstart_chim-1,iy+imstart_chim-1,:) & & + wrf_v(ix,iy) *wrf_uphi(ix+izstart_chim-1,iy+imstart_chim-1,:) enddo enddo do ime=1,nmerid do izo=1,nzonal np=domains%dom(rank+1)%izstart+izo+(domains%dom(rank+1)%imstart+ime-2)*nzonal_entire_domain-1 index_zonal=ixc(np)-izstart_chim+1 index_merid=iyc(np)-imstart_chim+1 if ((index_zonal.ne.izend_chim-izstart_chim+1).and.(index_merid.ne.imend_chim-imstart_chim+1) .and. (.not.iwrfgrid)) then w3dout(izo,ime,:) = (1.-wxc(np))*(1.-wyc(np))*w3d(index_zonal,index_merid,:) & + wxc(np) *(1.-wyc(np))*w3d(index_zonal+1,index_merid,:) & + (1.-wxc(np))* wyc(np) *w3d(index_zonal, index_merid+1,:) & + wxc(np) * wyc(np) *w3d(index_zonal+1,index_merid+1,:) else w3dout(izo,ime,:) = w3d(index_zonal,index_merid,:) endif chim_u(izo,ime,ks)=scalar(w3dout(izo,ime,:), & & chim_ulambda(izo+domains%dom(rank+1)%izstart-1,ime+domains%dom(rank+1)%imstart-1,:)) chim_v(izo,ime,ks)=scalar(w3dout(izo,ime,:), & & chim_uphi(izo+domains%dom(rank+1)%izstart-1,ime+domains%dom(rank+1)%imstart-1,:)) enddo enddo end subroutine wrf2chim_vec_2d !----------------------------------------------------------------------------------------------- subroutine wrf2chim_vec_3d(wrf_u,wrf_v,chim_u,chim_v,ks) ! inout because we want to respect slot 1 when filling slot 2 use chimere_common implicit none integer,intent(in) :: ks real(kind=iprec),intent(in) :: wrf_u(izend_chim-izstart_chim+1,imend_chim-imstart_chim+1,nlev_wrf) real(kind=iprec),intent(in) :: wrf_v(izend_chim-izstart_chim+1,imend_chim-imstart_chim+1,nlev_wrf) real(kind=iprec),intent(inout) :: chim_u(nzonal,nmerid,nlev_wrf,2),chim_v(nzonal,nmerid,nlev_wrf,2) integer :: iz do iz=1,nlev_wrf call wrf2chim_vec_2d(wrf_u(:,:,iz),wrf_v(:,:,iz),chim_u(:,:,iz,:),chim_v(:,:,iz,:),ks) enddo end subroutine wrf2chim_vec_3d function test_wrf_grid(nc_wrf,nx_wrf,ny_wrf) ! Routine tests wether the WRF horizontal grid is to the CHIMERE grid. The ! grids are considered identical if: ! - Grids have the same dimension ! - The center lons and lats of every WRF cell falls in the corresponding ! CHIMERE cell center ! Returns True if (and only if) both criteria above are met use ncvar use chimere_params, only : nzonal_entire_domain,nmerid_entire_domain use gridvars, only : xlati_entire_domain, xlong_entire_domain use interp_tools, only : calculate_corners, findcell implicit none ! Return value logical :: test_wrf_grid ! Input parameters type(ncfile),intent(IN) :: nc_WRF integer,intent(IN) :: nx_wrf,ny_wrf ! Local variables logical :: found real(kind=iprec),dimension(:,:,:,:),allocatable :: chim_corners_3d integer :: ixx,iyy,i,j,lstix,lstiy if ( nx_wrf.ne.nzonal_entire_domain .or. ny_wrf.ne.nmerid_entire_domain )then test_wrf_grid=.false. return endif ! Build the corners_3d array for findcell if (.not.allocated(chim_corners_3d)) allocate(chim_corners_3d(nzonal_entire_domain,nmerid_entire_domain,4,3)) call calculate_corners(xlong_entire_domain,xlati_entire_domain,chim_corners_3d,nzonal_entire_domain,nmerid_entire_domain) ! is every cell center of the WRF domain in the right CHIMERE cell ? do i=1,nzonal_entire_domain do j=1,nmerid_entire_domain lstix=i ! This is the expected result and therefore the best guess lstiy=j ! This is the expected result and therefore the best guess call findcell(nc_WRF%lon(i,j), nc_WRF%lat(i,j), ixx, iyy, & & nzonal_entire_domain, nmerid_entire_domain, chim_corners_3d, & & lstix, lstiy, found) if ( .not.(found) .or. (ixx.ne.i) .or. (iyy.ne.j) ) then test_wrf_grid=.false. return endif enddo enddo test_wrf_grid=.true. deallocate(chim_corners_3d) end function test_wrf_grid