PROGRAM arps2wrf,90 ! !################################################################## !################################################################## !###### ###### !###### PROGRAM ARPS2WRF ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Converts ADAS analysis file and/or EXT2ARPS dumps in ARPS grid to ! the corresponding WRF initialization file (wrfinput_d01) and ! the lateral boundary file (wrfbdy_d01) ! ! This program reads in a history file produced by ADAS and/or EXT2ARPS ! in any ARPS history format (except HDF 4 format), interpolates ! variables to WRF grid specified in the namelist file - arps2wrf.input ! and writes out these variables in NetCDF/PHDF5 format as wrfinput_d01 ! and wrfbdy_d01 ! ! Users can specify two options for the surface characteristics, such as ! vegetation type, soil type, and greenness fraction etc. ! ! Option 1 sfcinitopt = "ARPS", this program reads surface variables from ! the file created by ARPSSFC. the data must be at the same ! grid as that in the ADAS/EXT2ARPS data file. ! ! o Variables read in are: "soiltyp", "vegtyp", "veg"; ! o Variable, "xice", "xland", will be determined after the above ! variables are interpolated to WRF grid. ! o This option cannot set "SHDMAX","SHDMIN","SNOALB","ALBBCK". ! Should study whether it has negative effect to WRF run??? ! ! Option 2 sfcinitopt = "WRFSI", this program reads those variables from ! the static file 'static.wrfsi' created by gridgen_model of ! WRFSI package. The data must be at the same grid as those ! specified for WRF model in the input file, arps2wrf.input. ! ! (1). Users only need to run "window_domain_rt.pl" to generate ! WRFSI static file. (Actually, only gridgen_model.exe component). ! (2). One static file works for all runs in the same domain. ! ! ! INPUT ! arps2wrf.input NAMELIST ! runname.bin000000 ADAS time-dependent variable file ! runname.bingrdbas ADAS base state file ! runname.sfcdata ARPSSFC outputs ! OR ! static.wrfsi WRFSI static file ! ! OR ! outputs from ext2arps for lateral boundary conditions. ! ! OUTPUT ! wrfinput_d01 ! wrfbdy_d01 ! wrfnamelist.input ! ! Libraries: ! libarps.a ! libadas.a ! NetCDF library ! ! Subroutine calls: ! dtaread ! Subroutines defined in file interplib.f90 and wrf_iolib.f90 ! ! Use module ! wrf_metadata Defined in module_wrf_metadata.f90 ! !----------------------------------------------------------------------- ! ! AUTHOR: Yunheng Wang, Dan Weber, Keith Brewster ! 07/08/2003 ! ! MODIFICATION HISTORY: ! ! 03/10/2004 Yunheng Wang ! Upgraded to support WRF version 1.4 which was downloaded from ! WRF tiger team website on Feb. 10, 2004. ! ! 07/28/2004 Yunheng Wang ! Upgraded WRF V1.4 support to V2.0 support. So the unofficial release ! WRF V1.4 would not be supported any more since ARPS5.1.3. ! ! 11/20/2004 Yunheng Wang ! Reorangized to support PHDF5 format using WRF external IO-PHDF5 ! package. Upgraded to support WRFV2.0.3. Removed supports for ! previous WRF version. ! ! 01/10/2005 Yunheng Wang ! Added WRF internal binary support. Please note that WRF binary does ! not support random access. So the dumping order is important for ! binary data. The current code support WRFV2.0.3.1. The compatiblity ! between WRF versions will be bad, but it should not affect other ! IO format. ! !----------------------------------------------------------------------- ! ! ARPS DATA ARRAYS ! ! x x coordinate of grid points in physical/comp. space (m) ! y y coordinate of grid points in physical/comp. space (m) ! z z coordinate of grid points in computational space (m) ! zp z coordinate of grid points in physical space (m) ! zpsoil z coordinate of soil model (m) ! ! w vertical component of velocity in Cartesian ! coordinates (m/s). ! ! ptprt perturbation potential temperature (K) ! pprt perturbation pressure (Pascal) ! uprt perturbation x velocity component (m/s) ! vprt perturbation y velocity component (m/s) ! wprt perturbation z velocity component (m/s) ! ! qv water vapor mixing ratio (kg/kg) ! qc Cloud water mixing ratio (kg/kg) ! qr Rainwater mixing ratio (kg/kg) ! qi Cloud ice mixing ratio (kg/kg) ! qs Snow mixing ratio (kg/kg) ! qh Hail mixing ratio (kg/kg) ! ! ubar Base state x velocity component (m/s) ! vbar Base state y velocity component (m/s) ! wbar Base state z velocity component (m/s) ! ptbar Base state potential temperature (K) ! pbar Base state pressure (Pascal) ! rhobar Base state density (kg/m**3) ! qvbar Base state water vapor mixing ratio (kg/kg) ! ! soiltyp Soil type ! vegtyp Vegetation type ! lai Leaf Area Index ! roufns Surface roughness ! veg Vegetation fraction ! ! tsoil Soil Temperature (K) ! qsoil Soil moisture ! wetcanp Canopy water amount ! ! rain Total rain (raing + rainc) ! raing Grid supersaturation rain ! rainc Cumulus convective rain ! !----------------------------------------------------------------------- ! ! Variable Declarations: ! !----------------------------------------------------------------------- ! USE wrf_metadata ! WRF constants and metadata IMPLICIT NONE !---------------------------------------------------------------------- ! ! ARPS grid variables ! !--------------------------------------------------------------------- INTEGER :: nx,ny,nz ! Grid dimensions for ARPS. INTEGER :: nzsoil ! Soil levels INTEGER :: nstyps ! Maximum number of soil types. INTEGER, PARAMETER :: nhisfile_max = 100, nmax_domains = 100 INTEGER, PARAMETER :: max_vertical_levels = 100 CHARACTER(LEN=256) :: hisfile(nhisfile_max) CHARACTER(LEN=256) :: bdybasfn(nhisfile_max) INTEGER :: nhisfile,lengbf,lenfil CHARACTER(LEN=256) :: grdbasfn CHARACTER(LEN=256), DIMENSION(nmax_domains) :: adashisfn CHARACTER(LEN=256), DIMENSION(nmax_domains) :: adasbasfn INTEGER, DIMENSION(nmax_domains) :: hinfmt INTEGER :: max_dom LOGICAL, DIMENSION(nmax_domains) :: input_from_file INTEGER :: finexist(nhisfile_max) ! hisfile(1) and bdybasfn(1) are for WRF input file ! They can be any ARPS history dumps including ADAS output, ! ARPS output or EXT2ARPS output ! ! hisfile(2:nhisfil), bdybasfn(2:nhisfile) are for ! WRF lateral bounday files. They can be either ARPS history ! dumps or EXT2ARPS outputs CHARACTER(LEN=80) :: execname ! !----------------------------------------------------------------------- ! ! ARPS include files: ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' INCLUDE 'grid.inc' INCLUDE 'phycst.inc' INCLUDE 'mp.inc' ! !----------------------------------------------------------------------- ! ! Variables for mpi jobs ! !----------------------------------------------------------------------- INTEGER, PARAMETER :: fzone_arps = 3, fzone_wrf = 1 INTEGER :: ncompressx, ncompressy ! compression in x and y direction: ! ncompressx=nprocx_in/nproc_x ! ncompressy=nprocy_in/nproc_y INTEGER :: nxsm, nysm ! smaller domain INTEGER :: nxlg, nylg ! global domain INTEGER :: nxlg_wrf, nylg_wrf INTEGER :: nprocx_in, nprocy_in INTEGER :: ii,jj,ia,ja REAL, DIMENSION(:), POINTER :: xsm,ysm REAL, DIMENSION(:,:,:), POINTER :: zpsm,zpsoilsm REAL, DIMENSION(:,:,:), POINTER :: uprtsm, vprtsm, wprtsm, & ptprtsm, pprtsm, qvprtsm REAL, DIMENSION(:,:,:), POINTER :: qcsm, qrsm, qism, qssm, qhsm REAL, DIMENSION(:,:,:), POINTER :: ubarsm, vbarsm, wbarsm, & ptbarsm,pbarsm, qvbarsm,rhobarsm INTEGER, DIMENSION(:,:,:), POINTER :: soiltypsm INTEGER, DIMENSION(:,:), POINTER :: vegtypsm REAL, DIMENSION(:,:,:), POINTER :: stypfrctsm(:,:,:) REAL, DIMENSION(:,:,:,:), POINTER :: tsoilsm, qsoilsm REAL, DIMENSION(:,:,:), POINTER :: wetcanpsm REAL, DIMENSION(:,:), POINTER :: vegsm,snowdpthsm ! !----------------------------------------------------------------------- ! ! ARPS arrays to be read in: ! !----------------------------------------------------------------------- ! REAL, DIMENSION(:), POINTER :: x ! The x-coord. of the physical and ! computational grid. Defined at u-point. REAL, DIMENSION(:), POINTER :: y ! The y-coord. of the physical and ! computational grid. Defined at v-point. REAL, DIMENSION(:), POINTER :: z ! The z-coord. of the computational ! grid. Defined at w-point. REAL, DIMENSION(:,:,:), POINTER :: zp ! The height of the terrain. REAL, DIMENSION(:,:,:), POINTER :: zpsoil ! The height of the soil model. REAL, DIMENSION(:,:,:), POINTER :: uprt ! Perturbation u-velocity (m/s) REAL, DIMENSION(:,:,:), POINTER :: vprt ! Perturbation v-velocity (m/s) REAL, DIMENSION(:,:,:), POINTER :: wprt ! Perturbation w-velocity (m/s) REAL, DIMENSION(:,:,:), POINTER :: ptprt ! Perturbation potential temperature (K) REAL, DIMENSION(:,:,:), POINTER :: pprt ! Perturbation pressure (Pascal) REAL, DIMENSION(:,:,:), POINTER :: qvprt ! Perturbation water vapor specific ! humidity (kg/kg) REAL, DIMENSION(:,:,:), POINTER :: qc ! Cloud water mixing ratio (kg/kg) REAL, DIMENSION(:,:,:), POINTER :: qr ! Rain water mixing ratio (kg/kg) REAL, DIMENSION(:,:,:), POINTER :: qi ! Cloud ice mixing ratio (kg/kg) REAL, DIMENSION(:,:,:), POINTER :: qs ! Snow mixing ratio (kg/kg) REAL, DIMENSION(:,:,:), POINTER :: qh ! Hail mixing ratio (kg/kg) REAL, DIMENSION(:,:,:), POINTER :: qv ! Water vapor specific humidity (kg/kg) REAL, DIMENSION(:,:,:), POINTER :: ubar ! Base state u-velocity (m/s) REAL, DIMENSION(:,:,:), POINTER :: vbar ! Base state v-velocity (m/s) REAL, DIMENSION(:,:,:), POINTER :: wbar ! Base state w-velocity (m/s) REAL, DIMENSION(:,:,:), POINTER :: ptbar ! Base state potential temperature (K) REAL, DIMENSION(:,:,:), POINTER :: pbar ! Base state pressure (Pascal) REAL, DIMENSION(:,:,:), POINTER :: rhobar ! Base state air density (kg/m**3) REAL, DIMENSION(:,:,:), POINTER :: qvbar ! Base state water vapor specific INTEGER, DIMENSION(:,:,:), POINTER :: soiltyp ! Soil type REAL, DIMENSION(:,:,:), POINTER :: stypfrct ! Soil type INTEGER, DIMENSION(:,:), POINTER :: vegtyp ! Vegetation type REAL, DIMENSION(:,:), POINTER :: veg ! Vegetation fraction REAL, DIMENSION(:,:,:,:), POINTER :: tsoil ! Soil Temperature (K) REAL, DIMENSION(:,:,:,:), POINTER :: qsoil ! Soil Moisture REAL, DIMENSION(:,:,:), POINTER :: wetcanp ! Canopy water amount REAL, DIMENSION(:,:), POINTER :: snowdpth ! Snow depth (m) REAL, DIMENSION(:,:), POINTER :: dum2da REAL, DIMENSION(:,:,:), POINTER :: dum3da !----------------------------------------------------------------------- ! ! Namelist definitions for ARPS2WRF.input ! ! sfcdt Specify surface characteristics ! bdyspc Obtain boundary input files (ARPS format) ! wrf_grid Define WRF horizontal and vertical grid ! interp_options Choose interpolation scheme ! wrf_opts WRF options from namelist.input ! output Ouput options ! !----------------------------------------------------------------------- ! INTEGER, DIMENSION(nmax_domains) :: nxnm, nynm,nznm,nzsoilnm,nstypsnm INTEGER, DIMENSION(nmax_domains) :: nx_wrfnm, ny_wrfnm REAL, DIMENSION(nmax_domains) :: dx_wrfnm, dy_wrfnm INTEGER, DIMENSION(nmax_domains) :: i_parent_start, j_parent_start INTEGER, DIMENSION(nmax_domains) :: parent_id INTEGER :: nx_wrf ! = nx-2 if the same grid as ARPS INTEGER :: ny_wrf ! = ny-2 if the same grid as ARPS INTEGER :: nz_wrf ! = nz-2 if the same grid as ARPS ! All are staggered values INTEGER :: nzsoil_wrf ! WRF number of soil layers ! see sf_surface_physics REAL :: lattru_wrf(2) ! array of true latitude of WRF map projection REAL :: lontru_wrf ! true longitude of WRF map projection ! = trulon_wrf ! Namelist variable declaration CHARACTER(LEN=7), DIMENSION(nmax_domains) :: sfcinitopt ! either "ARPS" or "WRFSI" INTEGER :: create_bdy ! Create WRF boundary file INTEGER :: create_namelist ! Dump WRF namelist.input INTEGER, DIMENSION(nmax_domains) :: use_arps_grid ! Use ARPS horizontal grid as WRF grid INTEGER :: tintv_bdywrf ! The desired WRF boundary file interval INTEGER :: tintv_bdyin ! ARPS boundary file interval (in seconds) INTEGER :: mgrdbas ! Options for grid base file ! = 0 share same grid base as initial state file ! = 1 All ARPS boundary files share one grd base ! file but it is difference from the inital ! base file as specified using grdbasfn ! = 2 Each file has its own grid base file CHARACTER(LEN=256):: wrfnamelist ! file name for WRF namelist.input INTEGER :: mapproj_wrf ! Type of map projection in WRF model grid ! modproj = 1 Polar Stereographic ! projection. ! = 2 Mercator projection. ! = 3 Lambert projection. REAL :: sclfct_wrf ! Map scale factor. ! Distance on map, between two latitudes ! trulat1 and trulat2, ! is = (Distance on earth)*sclfct. ! For ARPS model runs, ! generally this is 1.0 REAL :: trulat1_wrf, trulat2_wrf, trulon_wrf ! 1st, 2nd real true latitude and true longitude ! of WRF map projection REAL :: ctrlat_wrf ! Center latitude of WRF model domain (deg. N) REAL :: ctrlon_wrf ! Center longitude of WRF model domain (deg. E) REAL :: dx_wrf ! WRF Grid spacing in x-direction REAL :: dy_wrf ! WRF Grid spacing in y-direction REAL, DIMENSION(nmax_domains) :: ptop ! WRF atmosphere top pressure in Pascal REAL :: zlevels_wrf(max_vertical_levels) ! WRF mass levels from 1.0 at surfact to ! 0.0 at atmosphere top INTEGER :: iorder ! order of polynomial for horizontal ! interpolation (1 or 2) INTEGER :: korder ! vertical interpolation order (1 or 2) INTEGER :: dyn_opt ! WRF dynamics option ! only works for = 2 Eulerian mass coordinate INTEGER :: diff_opt ! WRF diffusion option INTEGER :: km_opt ! WRF eddy coefficient option REAL, DIMENSION(nmax_domains) :: khdif ! Horizontal diffusion constant (m^2/s) REAL, DIMENSION(nmax_domains) :: kvdif ! Vertical diffusion constant (m^2/s) INTEGER,DIMENSION(nmax_domains) :: mp_physics ! WRF microphysics options != 2 Lin et al. scheme ! (QVAPOR,QRAIN,QSNOW,QCLOUD,QICE,QGRAUP) != 3 NCEP 3-class simple ice scheme ! (QVAPOR,QCLOUD,QICE,QRAIN,QSNOW) != 4 NCEP 5-class scheme ! (QVAPOR,QCLOUD,QICE,QRAIN,QSNOW) != 5 Ferrier (new Eta) microphysics ! (QVAPOR,QCLOUD) INTEGER, DIMENSION(nmax_domains) :: ra_lw_physics ! Longwave radiaiton option INTEGER, DIMENSION(nmax_domains) :: ra_sw_physics ! Shortwave radiation option INTEGER, DIMENSION(nmax_domains) :: sf_sfclay_physics ! WRF surface-layer option INTEGER, DIMENSION(nmax_domains) :: sf_surface_physics ! WRF land-surface option ! = 0 no land-surface ! (DO NOT use) ! = 1 Thermal diffusion scheme ! (nzsoil_wrf = 5) ! = 2 OSU land-surface model ! (nzsoil_wrf = 4) ! = 3 Do not use ! (nzsoil_wrf = 6) INTEGER, DIMENSION(nmax_domains) :: bl_pbl_physics ! boundary-layer option INTEGER, DIMENSION(nmax_domains) :: cu_physics ! cumulus option REAL :: dt ! time-step for advection INTEGER :: spec_bdy_width ! number of rows for specified boundary values nudging INTEGER :: nprocx_wrf ! Number of X direction processors for WRF run INTEGER :: nprocy_wrf ! Number of Y direction processors for WRF run INTEGER, DIMENSION(nmax_domains) :: frames_per_outfile INTEGER :: restart_interval REAL, DIMENSION(nmax_domains) :: radt,cudt INTEGER, DIMENSION(nmax_domains) :: parent_time_step_ratio INTEGER :: ifsnow, w_damping INTEGER :: io_form, io_form_history, io_form_restart INTEGER, DIMENSION(nmax_domains) :: history_interval LOGICAL :: pd_moist INTEGER :: qx_zero_out CHARACTER(LEN=256) :: staticdir, indir, outdir CHARACTER(LEN=4), PARAMETER :: fmtstr(7) = (/ 'bin_','xxxx','xxxx','xxxx','HDF5','xxxx','ncdf'/) ! !----------------------------------------------------------------------- ! ! ARPS working arrays ! !----------------------------------------------------------------------- REAL, ALLOCATABLE :: xs(:),ys(:) ! ARPS coordinate at scalar points REAL, ALLOCATABLE :: zps(:,:,:) ! ARPS vertical coordinate at scalar points REAL, ALLOCATABLE :: temp(:,:,:) ! ARPS temperature REAL, ALLOCATABLE :: tem1(:,:,:), tem2(:,:,:), tem3(:,:,:) !----------------------------------------------------------------------- ! ! WRF grid related variables ! !----------------------------------------------------------------------- REAL, ALLOCATABLE :: lat_wrf(:,:,:) ! WRF grid point latitudes REAL, ALLOCATABLE :: lon_wrf(:,:,:) ! WRF grid point lontitudes REAL, ALLOCATABLE :: msft_wrf(:,:) ! Map scale factor on mass grid REAL, ALLOCATABLE :: msfu_wrf(:,:) ! Map scale factor on u grid REAL, ALLOCATABLE :: msfv_wrf(:,:) ! Map scale factor on v grid REAL, ALLOCATABLE :: zlevels_half(:) ! WRF vertical mass points at half levels REAL, ALLOCATABLE :: zs_wrf(:) ! Depth of center of soil layers REAL, ALLOCATABLE :: dzs_wrf(:) ! Thickness of soil layers !----------------------------------------------------------------------- ! ! Horizontally interpolation arrays ! !----------------------------------------------------------------------- REAL, ALLOCATABLE :: x2d(:,:,:) REAL, ALLOCATABLE :: y2d(:,:,:) ! WRF coordinate at ARPS grid INTEGER, ALLOCATABLE :: iloc(:,:,:) ! i index of WRF point in ARPS array. INTEGER, ALLOCATABLE :: jloc(:,:,:) ! j index of WRF point in ARPS array. REAL, ALLOCATABLE :: dxfld(:,:) ! interpolation arrays calculated REAL, ALLOCATABLE :: rdxfld(:,:) ! in advance to speed up the REAL, ALLOCATABLE :: dyfld(:,:) ! interpolation REAL, ALLOCATABLE :: rdyfld(:,:) ! !----------------------------------------------------------------------- ! ! WRF work arrays ! (These arrays defined at WRF grid horizontally and ! vertically at ARPS physical heights) ! !----------------------------------------------------------------------- REAL, ALLOCATABLE :: t_tmp(:,:,:) REAL, ALLOCATABLE :: p_tmp(:,:,:) REAL, ALLOCATABLE :: qv_tmp(:,:,:) REAL, ALLOCATABLE :: zps_tmp(:,:,:) REAL, ALLOCATABLE :: etap(:,:,:,:) ! air mass at ARPS physical height vertically REAL, ALLOCATABLE :: tsoil_tmp(:,:,:) ! Defined at ARPS soil layers REAL, ALLOCATABLE :: qsoil_tmp(:,:,:) REAL, ALLOCATABLE :: zpsoil_tmp(:,:,:) REAL, ALLOCATABLE :: work3d(:,:,:) REAL, ALLOCATABLE :: work2d(:,:) !----------------------------------------------------------------------- ! ! WRF ARRAYS ! (Those are arrays defined at WRF grid both horizontally ! and vetically) ! !----------------------------------------------------------------------- REAL, ALLOCATABLE :: mu(:,:) ! WRF air mass at surface REAL, ALLOCATABLE :: hterain_wrf(:,:) ! WRF topograph REAL, ALLOCATABLE :: u_wrf(:,:,:) ! u-velocity (m/s) REAL, ALLOCATABLE :: v_wrf(:,:,:) ! v-velocity (m/s) REAL, ALLOCATABLE :: w_wrf(:,:,:) ! w-velocity (m/s) REAL, ALLOCATABLE :: ph_wrf(:,:,:) ! perturbation geopotential REAL, ALLOCATABLE :: phb_wrf(:,:,:) ! base state geopotential REAL, ALLOCATABLE :: pot_wrf(:,:,:) ! total potential. temp. REAL, ALLOCATABLE :: pt_wrf(:,:,:) ! perturbation potential. temp. REAL, ALLOCATABLE :: pt_init_wrf(:,:,:) REAL, ALLOCATABLE :: qv_wrf(:,:,:) ! Water vapor mixing ratio REAL, ALLOCATABLE :: p_wrf(:,:,:) ! perturbation pressure REAL, ALLOCATABLE :: pb_wrf(:,:,:) ! reference pressure REAL, ALLOCATABLE :: qc_wrf(:,:,:) REAL, ALLOCATABLE :: qr_wrf(:,:,:) REAL, ALLOCATABLE :: qi_wrf(:,:,:) REAL, ALLOCATABLE :: qs_wrf(:,:,:) REAL, ALLOCATABLE :: qg_wrf(:,:,:) REAL, ALLOCATABLE :: mup_wrf(:,:) ! perturbation air mass REAL, ALLOCATABLE :: mub_wrf(:,:) ! base air mass REAL, ALLOCATABLE :: tem1_wrf(:,:,:) ! WRF temporary arrays REAL, ALLOCATABLE :: tem2_wrf(:,:,:) REAL, ALLOCATABLE :: tem3_wrf(:,:,:) REAL, ALLOCATABLE :: tem4_wrf(:,:,:) INTEGER, ALLOCATABLE :: soiltyp_wrf(:,:),vegtyp_wrf(:,:) REAL, ALLOCATABLE :: vegfrct_wrf(:,:) REAL, ALLOCATABLE :: xice(:,:),xland(:,:) ! flags why they are reals? REAL, ALLOCATABLE :: sst_wrf(:,:) REAL, ALLOCATABLE :: hgt_wrf(:,:),tmn_wrf(:,:) REAL, ALLOCATABLE :: shdmin(:,:),shdmax(:,:),albbck(:,:),snoalb(:,:) REAL, ALLOCATABLE :: tslb_wrf(:,:,:) REAL, ALLOCATABLE :: smois_wrf(:,:,:) REAL, ALLOCATABLE :: snowh_wrf(:,:) REAL, ALLOCATABLE :: canwat_wrf(:,:) !---------------------------------------------------------------------- ! ! OUTPUT work arrays ! !---------------------------------------------------------------------- TYPE(wrf_global_metadata) :: global_meta REAL, ALLOCATABLE :: uatv_wrf(:,:,:) REAL, ALLOCATABLE :: vatu_wrf(:,:,:) REAL :: lat_ll(4), lat_lr(4), lat_ul(4), lat_ur(4) REAL :: lon_ll(4), lon_lr(4), lon_ul(4), lon_ur(4) !----------------------------------------------------------------------- ! ! Boundary arrays ! !----------------------------------------------------------------------- REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ubdy3dtemp1, vbdy3dtemp1, & tbdy3dtemp1, pbdy3dtemp1, qbdy3dtemp1 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET :: mbdy2dtemp1 REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: ubdy3dtemp2, vbdy3dtemp2, & tbdy3dtemp2,pbdy3dtemp2, qbdy3dtemp2 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET :: mbdy2dtemp2 REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: bdy3d1temp1,bdy3d1temp2, & bdy3d1temp3,bdy3d1temp4,bdy3d1temp5 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET :: bdy2d1temp1 REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET :: bdy3d2temp1,bdy3d2temp2, & bdy3d2temp3,bdy3d2temp4,bdy3d2temp5 REAL, DIMENSION(:,:), ALLOCATABLE, TARGET :: bdy2d2temp1 REAL, DIMENSION(:,:,:), POINTER :: ubdy3dtempc,vbdy3dtempc, & tbdy3dtempc,pbdy3dtempc,qbdy3dtempc REAL, DIMENSION(:,:), POINTER :: mbdy2dtempc REAL, DIMENSION(:,:,:), POINTER :: ubdy3dtempcn,vbdy3dtempcn, & tbdy3dtempcn,pbdy3dtempcn,qbdy3dtempcn REAL, DIMENSION(:,:), POINTER :: mbdy2dtempcn REAL, DIMENSION(:,:,:), ALLOCATABLE :: bdys,bdyn,bdyw,bdye REAL, DIMENSION(:,:,:), ALLOCATABLE :: blgs,blgn,blgw,blge !----------------------------------------------------------------------- ! ! Misc local variables ! !----------------------------------------------------------------------- ! INTEGER :: i,j,k,n,ifile INTEGER :: istatus, lenstr, ireturn INTEGER :: abstime, abstime1 INTEGER :: abstimep, abstimec, abstimen, abstimecn REAL :: time, gmthr INTEGER :: domid INTEGER :: abstdiff1, abstdiff2, abstdiff LOGICAL :: hinterp_needed LOGICAL :: fexist INTEGER :: count_bdy,count_in_loop REAL :: w1, w2 REAL :: latnot(2),ctrx, ctry, scalef REAL :: swx, swy REAL :: swx_wrf, swy_wrf CHARACTER(LEN=256), DIMENSION(nmax_domains) :: sfcdtfns CHARACTER(LEN=256), DIMENSION(nmax_domains) :: cdummy INTEGER, DIMENSION(nmax_domains) :: wrftrnopt CHARACTER(LEN=256) :: filename CHARACTER(LEN=24) :: times_str CHARACTER(LEN=24) :: nextbdytimes INTEGER :: nchr, ncid, nchout REAL :: rdummy ! INTEGER :: i_parent_end, j_parent_end INTEGER :: grid_ratio INTEGER, DIMENSION(nmax_domains) :: parent_grid_ratio REAL, DIMENSION(nmax_domains) :: xsub0, ysub0 INTEGER, DIMENSION(nmax_domains) :: year1,month1,day1,hour1,minute1,second1 INTEGER :: year2,month2,day2,hour2,minute2,second2 !WDT 2004-01-20 GMB: set SST ! REAL :: cnt,sst_sum,dist_weight ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! CALL mpinit_proc ! !----------------------------------------------------------------------- ! ! Initializations ! !----------------------------------------------------------------------- ! IF(myproc == 0) WRITE(6,'(14(/5x,a),/)') & '###############################################################',& '###############################################################',& '#### ####',& '#### Welcome to ARPS2WRF ####',& '#### ####',& '#### Program that reads in history files generated by ####',& '#### ADAS/ARPS and produces files for WRF to start: ####',& '#### ####',& '#### wrfinput_d01, wrfbdy_d01, namelist.input ####',& '#### ####',& '###############################################################',& '###############################################################' !----------------------------------------------------------------------- ! ! First, read in namelist input ! !----------------------------------------------------------------------- dyn_opt = 2 CALL readnamelist(0,max_dom,input_from_file, & hinfmt,adasbasfn,adashisfn,bdybasfn,hisfile,nhisfile,finexist, & nxnm,nynm,nznm,nzsoilnm,nstypsnm,nprocx_in,nprocy_in, & ncompressx,ncompressy, & use_arps_grid,nx_wrfnm,ny_wrfnm,nz_wrf,zlevels_wrf,ptop, & i_parent_start,j_parent_start,parent_id, & mapproj_wrf,sclfct_wrf,lattru_wrf,lontru_wrf, & ctrlat_wrf,ctrlon_wrf,dx_wrfnm,dy_wrfnm,dt, & base_pres,base_temp,base_lapse, & sfcinitopt,wrftrnopt,sfcdtfns,cdummy,cdummy,rdummy,rdummy, & create_bdy,mgrdbas,tintv_bdywrf,tintv_bdyin,spec_bdy_width, & diff_opt,km_opt,khdif,kvdif,mp_physics,ra_lw_physics, & ra_sw_physics,sf_sfclay_physics,sf_surface_physics,bl_pbl_physics, & cu_physics, nprocx_wrf,nprocy_wrf,frames_per_outfile, & restart_interval,radt,cudt,ifsnow,w_damping,parent_time_step_ratio,& iorder,korder,io_form,qx_zero_out,create_namelist,wrfnamelist, & io_form_history, io_form_restart,history_interval, & indir,outdir,staticdir,pd_moist,istatus) filename = ' ' execname = ' ' IF (mp_opt == 0) THEN execname = 'ARPS2WRF (serial)' ELSE WRITE(execname,'(a,2(I4.4,a))') 'ARPS2WRF (Parallel, ',nproc_x,'x',nproc_y,')' END IF !----------------------------------------------------------------------- ! ! Begin to loop for each domain ! !----------------------------------------------------------------------- DO domid = 1,max_dom IF (input_from_file(domid) .OR. domid == 1) THEN nx = nxnm(domid) ny = nynm(domid) nz = nznm(domid) nzsoil = nzsoilnm(domid) nstyps = nstypsnm(domid) nx_wrf = nx_wrfnm(domid) ny_wrf = ny_wrfnm(domid) dx_wrf = dx_wrfnm(domid) dy_wrf = dy_wrfnm(domid) IF (domid > 1) THEN create_bdy = 0 nhisfile = 1 hisfile (1) = adashisfn(domid) bdybasfn(1) = adasbasfn(domid) END IF WRITE(sfcdtfl,'(a)') sfcdtfns(domid) !----------------------------------------------------------------------- ! ! Setting microphysics variables ! !----------------------------------------------------------------------- CALL set_mp_physics_variables(domid,mp_physics(domid),istatus) IF (istatus /= 0) CALL arpsstop('ERROR in set_mp_physics_variables.',1) ! !----------------------------------------------------------------------- ! ! Allocate ARPS arrays ! !----------------------------------------------------------------------- !write(0,*) 'Deallocating' IF (ASSOCIATED(x)) THEN DEALLOCATE(x,y,z,zp,zpsoil) DEALLOCATE(uprt,vprt,wprt,ptprt,pprt,qvprt) DEALLOCATE(qv,qc,qr,qi,qs,qh) DEALLOCATE(ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar) DEALLOCATE(soiltyp,stypfrct,vegtyp,veg) DEALLOCATE(tsoil,qsoil,wetcanp,snowdpth) DEALLOCATE(dum2da, dum3da) END IF !write(0,*) 'reallocateing' ALLOCATE(x (nx)) ALLOCATE(y (ny)) ALLOCATE(z (nz)) ALLOCATE(zp (nx,ny,nz)) ALLOCATE(zpsoil (nx,ny,nzsoil)) ALLOCATE(uprt (nx,ny,nz)) ALLOCATE(vprt (nx,ny,nz)) ALLOCATE(wprt (nx,ny,nz)) ALLOCATE(ptprt (nx,ny,nz)) ALLOCATE(pprt (nx,ny,nz)) ALLOCATE(qvprt (nx,ny,nz)) ALLOCATE(qv (nx,ny,nz)) ALLOCATE(qc (nx,ny,nz)) ALLOCATE(qr (nx,ny,nz)) ALLOCATE(qi (nx,ny,nz)) ALLOCATE(qs (nx,ny,nz)) ALLOCATE(qh (nx,ny,nz)) ALLOCATE(ubar (nx,ny,nz)) ALLOCATE(vbar (nx,ny,nz)) ALLOCATE(wbar (nx,ny,nz)) ALLOCATE(ptbar (nx,ny,nz)) ALLOCATE(pbar (nx,ny,nz)) ALLOCATE(rhobar (nx,ny,nz)) ALLOCATE(qvbar (nx,ny,nz)) ALLOCATE(soiltyp (nx,ny,nstyps)) ALLOCATE(stypfrct(nx,ny,nstyps)) ALLOCATE(vegtyp (nx,ny)) ALLOCATE(veg (nx,ny)) ALLOCATE(tsoil (nx,ny,nzsoil,0:nstyps)) ALLOCATE(qsoil (nx,ny,nzsoil,0:nstyps)) ALLOCATE(wetcanp(nx,ny,0:nstyps)) ALLOCATE(snowdpth(nx,ny)) ALLOCATE(dum2da(nx,ny)) ALLOCATE(dum3da(nx,ny,nz)) x =0.0 y =0.0 z =0.0 zp =0.0 zpsoil =0.0 uprt =0.0 vprt =0.0 wprt =0.0 ptprt =0.0 pprt =0.0 qvprt =0.0 qc =0.0 qr =0.0 qi =0.0 qs =0.0 qh =0.0 ubar =0.0 vbar =0.0 wbar =0.0 ptbar =0.0 pbar =0.0 rhobar =0.0 qvbar =0.0 qv =0.0 soiltyp =0.0 stypfrct=0.0 vegtyp =0.0 veg =0.0 tsoil =0.0 qsoil =0.0 wetcanp =0.0 snowdpth=0.0 dum2da = 0.0 !write(0,*) 'second reallocating' IF (ALLOCATED(xs)) DEALLOCATE(xs,ys,zps,temp) ALLOCATE(xs(nx), STAT = istatus) ALLOCATE(ys(ny), STAT = istatus) ALLOCATE(zps(nx,ny,nz), STAT = istatus) ALLOCATE(temp(nx,ny,nz), STAT = istatus) xs = 0.0 ys = 0.0 zps = 0.0 temp = 0.0 !--------------------------------------------------------------------------- ! ! Allocate WRF arrays ! !--------------------------------------------------------------------------- IF (ALLOCATED(zlevels_half)) THEN !write(0,*) 'deallocating zlevels_half' DEALLOCATE(zlevels_half) !write(0,*) 'deallocating mu' DEALLOCATE(mu,hterain_wrf) !write(0,*) 'deallocating integer' DEALLOCATE(soiltyp_wrf,vegtyp_wrf) !write(0,*) 'deallocating next' DEALLOCATE(vegfrct_wrf,xice,xland) !write(0,*) 'deallocating sst_wrf' DEALLOCATE(sst_wrf,hgt_wrf,tmn_wrf,shdmin,shdmax,albbck,snoalb) !write(0,*) 'deallocating snowh_wrf' DEALLOCATE(snowh_wrf,canwat_wrf) END IF !write(0,*) 'allocating wrf arrys' ALLOCATE(zlevels_half(nz_wrf-1), STAT= istatus) DO k = 1, nz_wrf-1 zlevels_half(k) = (zlevels_wrf(k)+zlevels_wrf(k+1))*0.5 END DO ALLOCATE(mu (nx_wrf,ny_wrf), STAT = istatus) ALLOCATE(hterain_wrf(nx_wrf,ny_wrf), STAT = istatus) ALLOCATE(soiltyp_wrf(nx_wrf,ny_wrf), STAT = istatus) !write(0,*) 'allocate status of soiltyp_wrf ',istatus ALLOCATE(vegtyp_wrf (nx_wrf,ny_wrf), STAT = istatus) ALLOCATE(vegfrct_wrf(nx_wrf,ny_wrf), STAT = istatus) ALLOCATE(xice (nx_wrf,ny_wrf), STAT = istatus) ALLOCATE(xland (nx_wrf,ny_wrf), STAT = istatus) ALLOCATE(sst_wrf(nx_wrf, ny_wrf ), STAT = istatus) ALLOCATE(hgt_wrf(nx_wrf, ny_wrf ), STAT = istatus) ALLOCATE(tmn_wrf(nx_wrf, ny_wrf ), STAT = istatus) ALLOCATE(shdmin (nx_wrf, ny_wrf ), STAT = istatus) ALLOCATE(shdmax (nx_wrf, ny_wrf ), STAT = istatus) ALLOCATE(albbck (nx_wrf, ny_wrf ), STAT = istatus) ALLOCATE(snoalb (nx_wrf, ny_wrf ), STAT = istatus) ALLOCATE(snowh_wrf (nx_wrf,ny_wrf), STAT = istatus) ALLOCATE(canwat_wrf(nx_wrf,ny_wrf), STAT = istatus) IF (ALLOCATED(u_wrf)) THEN DEALLOCATE(u_wrf,v_wrf,w_wrf,ph_wrf,phb_wrf,pot_wrf) DEALLOCATE(pt_wrf,p_wrf,pb_wrf,qv_wrf,qc_wrf,pt_init_wrf) DEALLOCATE(qr_wrf,qi_wrf,qs_wrf,qg_wrf,mup_wrf,mub_wrf) DEALLOCATE(tem1_wrf,tem2_wrf,tem3_wrf,tem4_wrf) DEALLOCATE(zs_wrf,dzs_wrf,msft_wrf,msfu_wrf,msfv_wrf) END IF ALLOCATE(u_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(v_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(w_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(ph_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(phb_wrf(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(pot_wrf(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(pt_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(pt_init_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(p_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(pb_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(qv_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(qc_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(qr_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(qi_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(qs_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(qg_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(mup_wrf(nx_wrf,ny_wrf), STAT = istatus) ALLOCATE(mub_wrf(nx_wrf,ny_wrf), STAT = istatus) ALLOCATE(tem1_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(tem2_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(tem3_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(tem4_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(zs_wrf (9), STAT = istatus) ! maximum soil layer is 9 ALLOCATE(dzs_wrf(9), STAT = istatus) ALLOCATE(msft_wrf(nx_wrf,ny_wrf), STAT = istatus) ALLOCATE(msfu_wrf(nx_wrf,ny_wrf), STAT = istatus) ALLOCATE(msfv_wrf(nx_wrf,ny_wrf), STAT = istatus) IF (ALLOCATED(lat_wrf)) DEALLOCATE(lat_wrf,lon_wrf) ALLOCATE(lat_wrf(nx_wrf,ny_wrf,3), STAT = istatus) ALLOCATE(lon_wrf(nx_wrf,ny_wrf,3), STAT = istatus) mu = 0.0 hterain_wrf = 0.0 soiltyp_wrf = 0 vegtyp_wrf = 0 xice = 0.0 xland = 0.0 u_wrf = 0.0 v_wrf = 0.0 ph_wrf = 0.0 phb_wrf = 0.0 pot_wrf = 0.0 pt_wrf = 0.0 pt_init_wrf = 0.0 p_wrf = 0.0 pb_wrf = 0.0 qv_wrf = 0.0 qc_wrf = 0.0 qr_wrf = 0.0 qi_wrf = 0.0 qs_wrf = 0.0 qg_wrf = 0.0 mup_wrf = 0.0 mub_wrf = 0.0 zs_wrf = 0.0 dzs_wrf = 0.0 snowh_wrf = 0.0 canwat_wrf = 0.0 !----------------------------------------------------------------------- ! ! Allocate temporary arrays ! !----------------------------------------------------------------------- nxlg_wrf = (nx_wrf-fzone_wrf)*nproc_x+fzone_wrf nylg_wrf = (ny_wrf-fzone_wrf)*nproc_y+fzone_wrf IF (ALLOCATED(tem1)) THEN DEALLOCATE(tem1,tem2,tem3) DEALLOCATE(t_tmp,p_tmp,qv_tmp,zps_tmp) DEALLOCATE(etap,work2d,work3d) END IF !write(0,*) 'allocating temporary arrys' ALLOCATE(tem1(nx,ny,nz), STAT= istatus) ALLOCATE(tem2(nx,ny,nz), STAT= istatus) ALLOCATE(tem3(nx,ny,nz), STAT= istatus) ALLOCATE(t_tmp (nx_wrf,ny_wrf,nz), STAT = istatus) ALLOCATE(p_tmp (nx_wrf,ny_wrf,nz), STAT = istatus) ALLOCATE(qv_tmp (nx_wrf,ny_wrf,nz), STAT = istatus) ALLOCATE(zps_tmp (nx_wrf,ny_wrf,nz), STAT = istatus) ALLOCATE(etap (nx_wrf,ny_wrf,nz,3), STAT = istatus) ALLOCATE(work3d(nx_wrf,ny_wrf,nz), STAT = istatus) ALLOCATE(work2d(nx_wrf,ny_wrf), STAT = istatus) !write(0,*) 'create_bdy = ',create_bdy IF(create_bdy >= 1) THEN ALLOCATE(ubdy3dtemp1(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(vbdy3dtemp1(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(tbdy3dtemp1(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(pbdy3dtemp1(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(qbdy3dtemp1(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(mbdy2dtemp1(nx_wrf,ny_wrf), STAT = istatus) ALLOCATE(ubdy3dtemp2(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(vbdy3dtemp2(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(tbdy3dtemp2(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(pbdy3dtemp2(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(qbdy3dtemp2(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(mbdy2dtemp2(nx_wrf,ny_wrf), STAT = istatus) ALLOCATE(bdys(nx_wrf,nz_wrf, spec_bdy_width), STAT = istatus) ALLOCATE(bdyn(nx_wrf,nz_wrf, spec_bdy_width), STAT = istatus) ALLOCATE(bdyw(ny_wrf,nz_wrf, spec_bdy_width), STAT = istatus) ALLOCATE(bdye(ny_wrf,nz_wrf, spec_bdy_width), STAT = istatus) ALLOCATE(blgs(nxlg_wrf, nz_wrf, spec_bdy_width), STAT = istatus) ALLOCATE(blgn(nxlg_wrf, nz_wrf, spec_bdy_width), STAT = istatus) ALLOCATE(blgw(nylg_wrf, nz_wrf, spec_bdy_width), STAT = istatus) ALLOCATE(blge(nylg_wrf, nz_wrf, spec_bdy_width), STAT = istatus) blgs = 0.0 blgn = 0.0 blgw = 0.0 blge = 0.0 bdys = 0.0 bdyn = 0.0 bdyw = 0.0 bdye = 0.0 END IF IF (ALLOCATED(x2d)) DEALLOCATE(x2d,y2d,iloc,jloc) ALLOCATE(x2d (nx_wrf,ny_wrf,3), STAT = istatus) ALLOCATE(y2d (nx_wrf,ny_wrf,3), STAT = istatus) ALLOCATE(iloc (nx_wrf,ny_wrf,3), STAT = istatus) ALLOCATE(jloc (nx_wrf,ny_wrf,3), STAT = istatus) x2d = 0.0 y2d = 0.0 iloc = 0 jloc = 0 ! ! mpi variables ! nxsm = (nx-fzone_arps)/ncompressx + fzone_arps nysm = (ny-fzone_arps)/ncompressy + fzone_arps nxlg = (nx-fzone_arps)*nproc_x + fzone_arps nylg = (ny-fzone_arps)*nproc_y + fzone_arps IF(ncompressx > 1 .OR. ncompressy > 1) THEN IF (ASSOCIATED(xsm)) THEN DEALLOCATE(xsm,ysm,zpsm,zpsoilsm) DEALLOCATE(uprtsm,vprtsm,wprtsm,ptprtsm,pprtsm,qvprtsm) DEALLOCATE(qcsm,qrsm,qism,qssm,qhsm) DEALLOCATE(ubarsm,vbarsm,wbarsm,ptbarsm,pbarsm,rhobarsm,qvbarsm) DEALLOCATE(soiltypsm,stypfrctsm,vegtypsm,vegsm) DEALLOCATE(tsoilsm,qsoilsm,wetcanpsm,snowdpthsm) END IF ALLOCATE(xsm (nxsm),STAT=istatus) ALLOCATE(ysm (nysm),STAT=istatus) ALLOCATE(zpsm (nxsm,nysm,nz), STAT=istatus) ALLOCATE(zpsoilsm(nxsm,nysm,nzsoil),STAT=istatus) ALLOCATE(uprtsm(nxsm,nysm,nz), STAT=istatus) ALLOCATE(vprtsm(nxsm,nysm,nz), STAT=istatus) ALLOCATE(wprtsm(nxsm,nysm,nz), STAT=istatus) ALLOCATE(ptprtsm(nxsm,nysm,nz),STAT=istatus) ALLOCATE(pprtsm(nxsm,nysm,nz), STAT=istatus) ALLOCATE(qvprtsm(nxsm,nysm,nz),STAT=istatus) ALLOCATE(qcsm (nxsm,nysm,nz),STAT=istatus) ALLOCATE(qrsm (nxsm,nysm,nz),STAT=istatus) ALLOCATE(qism (nxsm,nysm,nz),STAT=istatus) ALLOCATE(qssm (nxsm,nysm,nz),STAT=istatus) ALLOCATE(qhsm (nxsm,nysm,nz),STAT=istatus) ALLOCATE(ubarsm (nxsm,nysm,nz),STAT=istatus) ALLOCATE(vbarsm (nxsm,nysm,nz),STAT=istatus) ALLOCATE(wbarsm (nxsm,nysm,nz),STAT=istatus) ALLOCATE(ptbarsm (nxsm,nysm,nz),STAT=istatus) ALLOCATE(pbarsm (nxsm,nysm,nz),STAT=istatus) ALLOCATE(rhobarsm(nxsm,nysm,nz),STAT=istatus) ALLOCATE(qvbarsm (nxsm,nysm,nz),STAT=istatus) ALLOCATE(soiltypsm (nxsm,nysm,nstyps),STAT=istatus) ALLOCATE(stypfrctsm(nxsm,nysm,nstyps),STAT=istatus) ALLOCATE(vegtypsm (nxsm,nysm),STAT=istatus) ALLOCATE(vegsm (nxsm,nysm),STAT=istatus) ALLOCATE(tsoilsm (nxsm,nysm,nzsoil,0:nstyps),STAT=istatus) ALLOCATE(qsoilsm (nxsm,nysm,nzsoil,0:nstyps),STAT=istatus) ALLOCATE(wetcanpsm (nxsm,nysm,0:nstyps),STAT=istatus) ALLOCATE(snowdpthsm(nxsm,nysm),STAT=istatus) CALL check_alloc_status(istatus, "arpsplt:snowdpth") xsm = 0.0 ysm = 0.0 zpsm = 0.0 zpsoilsm= 0.0 uprtsm = 0.0 vprtsm = 0.0 wprtsm = 0.0 ptprtsm = 0.0 pprtsm = 0.0 qvprtsm = 0.0 qcsm = 0.0 qrsm = 0.0 qism = 0.0 qssm = 0.0 qhsm = 0.0 ubarsm = 0.0 vbarsm = 0.0 wbarsm = 0.0 ptbarsm = 0.0 pbarsm = 0.0 rhobarsm= 0.0 qvbarsm = 0.0 soiltypsm = 0.0 stypfrctsm= 0.0 vegtypsm = 0.0 vegsm = 0.0 tsoilsm = 0.0 qsoilsm = 0.0 wetcanpsm = 0.0 snowdpthsm= 0.0 ELSE xsm => x ysm => y zpsm => zp zpsoilsm => zpsoil uprtsm => uprt vprtsm => vprt wprtsm => wprt ptprtsm => ptprt pprtsm => pprt qvprtsm => qvprt qcsm => qc qrsm => qr qism => qi qssm => qs qhsm => qh ubarsm => ubar vbarsm => vbar wbarsm => wbar ptbarsm => ptbar pbarsm => pbar rhobarsm => rhobar qvbarsm => qvbar soiltypsm => soiltyp stypfrctsm => stypfrct vegtypsm => vegtyp vegsm => veg tsoilsm => tsoil qsoilsm => qsoil wetcanpsm => wetcanp snowdpthsm => snowdpth END IF !-------------------------------------------------------------------- ! ! Begin loop over all of the history files ! !-------------------------------------------------------------------- count_bdy = 0 DO ifile = 1,nhisfile !write(0,*) 'domid,fileid ',domid,ifile lenfil = LEN_TRIM(hisfile(ifile)) CALL strlnth( hisfile(ifile), lenfil) IF (finexist(ifile) == 1) THEN ! file to be read exist DO jj = 1, ncompressy DO ii = 1, ncompressx IF (mp_opt > 0 .AND. readsplit <= 0) THEN WRITE(hisfile(ifile)(lenfil+1:lenfil+5),'(a,i2.2,i2.2)') '_', & ncompressx*(loc_x-1)+ii,ncompressy*(loc_y-1)+jj lenfil= lenfil + 5 WRITE(6,'(a,I5,2a/)') 'Processor: ',myproc, & ' reading file: ', hisfile(ifile)(1:lenfil) ELSE IF (myproc == 0) THEN WRITE(6,'(/2a/)') ' Reading file: ',hisfile(ifile)(1:lenfil) END IF ! !----------------------------------------------------------------------- ! ! Set the gridread parameter to 0 so that the boundary ! grid/base file will be read. ! !----------------------------------------------------------------------- ! IF( mgrdbas == 0 .OR. ifile == 1) THEN grdbasfn = bdybasfn(1) lengbf = LEN_TRIM(grdbasfn) IF (mp_opt > 0 .AND. readsplit <= 0) THEN WRITE(grdbasfn(lengbf+1:lengbf+5),'(a,i2.2,i2.2)') '_', & ncompressx*(loc_x-1)+ii,ncompressy*(loc_y-1)+jj lengbf= lengbf + 5 END IF IF (myproc == 0) WRITE(6,'(1x,a,a)') 'The grid/base name is ', grdbasfn(1:lengbf) IF (ifile == 1) CALL setgbrd(0) ELSE IF (mgrdbas == 1 .AND. ifile == 2) THEN grdbasfn = bdybasfn(2) lengbf = LEN_TRIM(grdbasfn) IF (mp_opt > 0 .AND. readsplit <= 0) THEN WRITE(grdbasfn(lengbf+1:lengbf+5),'(a,i2.2,i2.2)') '_', & ncompressx*(loc_x-1)+ii,ncompressy*(loc_y-1)+jj lengbf= lengbf + 5 END IF IF (myproc == 0) WRITE(6,'(1x,a,a)') ' The grid/base name is ', grdbasfn(1:lengbf) CALL setgbrd(0) ELSE IF (mgrdbas == 2 .AND. ifile >= 2) THEN grdbasfn = bdybasfn(ifile) lengbf = LEN_TRIM(grdbasfn) IF (mp_opt > 0 .AND. readsplit <= 0) THEN WRITE(grdbasfn(lengbf+1:lengbf+5),'(a,i2.2,i2.2)') '_', & ncompressx*(loc_x-1)+ii,ncompressy*(loc_y-1)+jj lengbf= lengbf + 5 END IF IF (myproc == 0) WRITE(6,'(1x,a,a)') ' The grid/base name is ', grdbasfn(1:lengbf) CALL setgbrd(0) END IF ! !----------------------------------------------------------------------- ! ! Read all input data arrays ! !----------------------------------------------------------------------- ! CALL dtaread(nx,ny,nz,nzsoil,nstyps,hinfmt(domid),nchr,grdbasfn(1:lengbf), & lengbf, hisfile(ifile)(1:lenfil),lenfil,time, & xsm,ysm,z,zpsm,zpsoilsm,uprtsm,vprtsm,wprtsm,ptprtsm,pprtsm, & qvprtsm, qcsm, qrsm, qism, qssm, qhsm, dum3da,dum3da,dum3da, & ubarsm, vbarsm, wbarsm, ptbarsm, pbarsm, rhobarsm, qvbarsm, & soiltypsm,stypfrctsm,vegtypsm,dum2da,dum2da,vegsm, & tsoilsm,qsoilsm,wetcanpsm,snowdpthsm, & dum2da,dum2da,dum3da, & dum3da,dum2da,dum2da,dum2da,dum2da, & dum2da,dum2da,dum2da,dum2da, & ireturn, tem1,tem2,tem3) !----------------------------------------------------------------------- ! ! ireturn = 0 for a successful read ! !----------------------------------------------------------------------- ! IF( ireturn /= 0 ) CALL arpsstop('dtaread errors.',1) IF(ncompressx > 1 .OR. ncompressy > 1) THEN ! need join DO j = 1, nysm ja = (jj-1)*(nysm-3)+j DO i = 1, nxsm ia = (ii-1)*(nxsm-3)+i x(ia) = xsm(i) vegtyp(ia,ja) = vegtypsm(i,j) veg(ia,ja) = vegsm(i,j) snowdpth(ia,ja) = snowdpthsm(i,j) zp(ia,ja,:) = zpsm(i,j,:) zpsoil(ia,ja,:) = zpsoilsm(i,j,:) uprt(ia,ja,:) = uprtsm(i,j,:) vprt(ia,ja,:) = vprtsm(i,j,:) wprt(ia,ja,:) = wprtsm(i,j,:) ptprt(ia,ja,:) = ptprtsm(i,j,:) pprt(ia,ja,:) = pprtsm(i,j,:) qvprt(ia,ja,:) = qvprtsm(i,j,:) qc(ia,ja,:) = qcsm(i,j,:) qr(ia,ja,:) = qrsm(i,j,:) qi(ia,ja,:) = qism(i,j,:) qs(ia,ja,:) = qssm(i,j,:) qh(ia,ja,:) = qhsm(i,j,:) ubar(ia,ja,:) = ubarsm(i,j,:) vbar(ia,ja,:) = vbarsm(i,j,:) wbar(ia,ja,:) = wbarsm(i,j,:) ptbar(ia,ja,:) = ptbarsm(i,j,:) pbar(ia,ja,:) = pbarsm(i,j,:) rhobar(ia,ja,:) = rhobarsm(i,j,:) qvbar(ia,ja,:) = qvbarsm(i,j,:) soiltyp(ia,ja,:) = soiltypsm(i,j,:) stypfrct(ia,ja,:) = stypfrctsm(i,j,:) tsoil(ia,ja,:,:) = tsoilsm(i,j,:,:) qsoil(ia,ja,:,:) = qsoilsm(i,j,:,:) wetcanp(ia,ja,:) = wetcanpsm(i,j,:) END DO !i y(ja) = ysm(j) END DO !j END IF ! join needed END DO END DO ELSE ! file to be read does not exist WRITE(6,'(1x,3a)') 'WARNING: File ',hisfile(ifile)(1:lenfil), & ' is skipped because it does not exist.' CYCLE ! read next file, time interpolation needed END IF ! Data read finished CALL ctim2abss( year,month,day,hour,minute,second, abstime ) abstime = abstime + INT(time) CALL abss2ctim( abstime, year, month, day, hour, minute, second ) CALL julday( year, month, day, jday) abstimep = abstimen abstimen = abstime !======================================================================= ! ! Only do the following for the first time level ! !======================================================================= IF( ifile == 1 ) THEN gmthr = hour + minute/60. + second/3600. IF(use_arps_grid(domid) == 1) THEN mapproj_wrf = mapproj sclfct_wrf = sclfct trulat1_wrf = trulat1 trulat2_wrf = trulat2 trulon_wrf = trulon ctrlat_wrf = ctrlat ctrlon_wrf = ctrlon dx_wrf = x(3)-x(2) dy_wrf = y(3)-y(2) hinterp_needed = .FALSE. dx_wrfnm(domid) = dx_wrf dy_wrfnm(domid) = dy_wrf ELSE trulat1_wrf = lattru_wrf(1) trulat2_wrf = lattru_wrf(2) trulon_wrf = lontru_wrf END IF IF (myproc == 0) THEN WRITE(6,'(/a/)') 'The ARPS grid and the WRF grid are:' WRITE(6,'(a38)') ' ARPS grid WRF grid' WRITE(6,'(a38)') ' ========== ==========' WRITE(6,'(a,I10,4x,I10)') ' nx = ',nx,nx_wrf WRITE(6,'(a,I10,4x,I10)') ' ny = ',ny,ny_wrf WRITE(6,'(a,I10,4x,I10)') ' mapproj = ',mapproj,mapproj_wrf WRITE(6,'(a,F10.2,4x,F10.2)') ' sclfct = ',sclfct,sclfct_wrf WRITE(6,'(a,F10.2,4x,F10.2)') ' trulat1 = ',trulat1,trulat1_wrf WRITE(6,'(a,F10.2,4x,F10.2)') ' trulat2 = ',trulat2,trulat2_wrf WRITE(6,'(a,F10.2,4x,F10.2)') ' trulon = ',trulon,trulon_wrf WRITE(6,'(a,F10.2,4x,F10.2)') ' ctrlat = ',ctrlat,ctrlat_wrf WRITE(6,'(a,F10.2,4x,F10.2)') ' ctrlon = ',ctrlon,ctrlon_wrf WRITE(6,'(a,F10.0,4x,F10.0)') ' dx = ',dx,dx_wrf WRITE(6,'(a,F10.0,4x,F10.0)') ' dy = ',dy,dy_wrf END IF IF(mapproj == mapproj_wrf .AND. & trulat1 == trulat1_wrf .AND. trulat2 == trulat2_wrf .AND. & trulon == trulon_wrf .AND. & ctrlat == ctrlat_wrf .AND. ctrlon == ctrlon_wrf .AND. & nx == nx_wrf + 2 .AND. ny == ny_wrf + 2 .AND. & ABS(dx_wrf-dx) < a_small_real_number .AND. ABS(dy_wrf-dy) < a_small_real_number ) THEN ! ARPS and WRF are at the same grid ! Horizontal interpolations are not necessary hinterp_needed = .FALSE. IF(myproc == 0) & WRITE(6,'(/a/)') ' NOTE: No horizontal interpolation will be performed.' ELSE hinterp_needed = .TRUE. WRITE(6,'(/a/)') ' NOTE: Horizontal interpolation will be performed.' END IF IF (hinterp_needed .AND. mp_opt > 0) THEN WRITE(6,'(1x,/2a/)') 'ARPS2WRF_mpi does not support horizontal ', & 'interpolation still. Please check and try again.' CALL arpsstop('mpi mode still does not implemented.',1) END IF ! ! It was put here becasue dx_wrf/dy_wrf is finalized only until now. ! CALL get_check_grid_ratio(domid,nx_wrf,ny_wrf,dx_wrf,dy_wrf, & dx_wrfnm(parent_id(domid)),dy_wrfnm(parent_id(domid)), & a_small_real_number,grid_ratio,istatus) IF (istatus /= 0) CALL arpsstop('Unacceptable parent_grid_ratio or nesting grid size.',1) parent_grid_ratio(domid) = grid_ratio year1(domid) = year ! Save the initial fields time month1(domid) = month day1(domid) = day hour1(domid) = hour minute1(domid)= minute second1(domid)= second year2 = year ! Just initialize the last boundary time month2 = month day2 = day hour2 = hour minute2 = minute second2 = second abstime1 = abstime abstimec = abstime ! !----------------------------------------------------------------------- ! ! Establish coordinate for ARPS scalar fields. ! !----------------------------------------------------------------------- ! DO i=1,nx-1 xs(i)=0.5*(x(i)+x(i+1)) END DO xs(nx)=2.*xs(nx-1)-xs(nx-2) DO j=1,ny-1 ys(j)=0.5*(y(j)+y(j+1)) END DO ys(ny)=2.*ys(ny-1)-ys(ny-2) DO k=1,nz-1 DO j=1,ny DO i=1,nx zps(i,j,k)=0.5*(zp(i,j,k)+zp(i,j,k+1)) END DO END DO END DO zps(:,:,nz) = 2*zps(:,:,nz-1)-zps(:,:,nz-2) ! !----------------------------------------------------------------------- ! ! Establish WRF map projection and find latitude and longitude of WRF grid. ! and compute the lat/lon at domain corners ! ! 1 -- T point, 2 -- U point, 3 -- V point, 4 -- massless point ! !----------------------------------------------------------------------- ! lattru_wrf(1) = trulat1_wrf lattru_wrf(2) = trulat2_wrf lontru_wrf = trulon_wrf ! sclfct_wrf changes inside CALL build_wrf_grid(domid,parent_id(domid),use_arps_grid, & dx_wrfnm(parent_id(domid)),dy_wrfnm(parent_id(domid)), & xsub0(parent_id(domid)),ysub0(parent_id(domid)),& i_parent_start(domid),j_parent_start(domid), & nx_wrf,ny_wrf,nxlg_wrf,nylg_wrf,dx_wrf,dy_wrf, & mapproj_wrf,sclfct_wrf,lattru_wrf,lontru_wrf, & ctrlat_wrf,ctrlon_wrf,swx_wrf,swy_wrf, & xsub0(domid),ysub0(domid), & lat_wrf,lon_wrf,lat_ll,lat_ul,lat_ur,lat_lr, & lon_ll,lon_ul,lon_ur,lon_lr,istatus) ! i_parent_end = i_parent_start(domid) + (nx_wrf-1)/grid_ratio ! j_parent_end = j_parent_start(domid) + (ny_wrf-1)/grid_ratio ! ! IF (i_parent_end > nx_wrfnm(parent_id(domid)) .OR. & ! j_parent_end > ny_wrfnm(parent_id(domid)) ) THEN ! WRITE(6,'(1x,2(a,I2),a)') 'Subdomain ',domid,' exceed its parent domain ', & ! parent_id(domid),' in size.' ! WRITE(6,'(1x,4(a,I8),a)') 'i_parent_start-end, j_parent_start-end = ', & ! i_parent_start(domid),' - ',i_parent_end,', ', & ! j_parent_start(domid),' - ',j_parent_end,'.' ! WRITE(6,'(1x,a,I2,a,2(I8,a),/)') 'Parent domain ',parent_id(domid),' size: ', & ! nx_wrfnm(parent_id(domid)),', ', & ! ny_wrfnm(parent_id(domid)),'.' ! ! CALL arpsstop('Subdomain is too large.',1) ! END IF !----------------------------------------------------------------------- ! ! Find x,y locations of WRF grid in terms of the ARPS grid. ! !----------------------------------------------------------------------- IF (hinterp_needed .OR. sfcinitopt(domid) == 'ARPS') THEN latnot(1) = trulat1 latnot(2) = trulat2 scalef = 1.0/sclfct IF (ALLOCATED(dxfld)) DEALLOCATE(dxfld,dyfld,rdxfld,rdyfld) ALLOCATE(dxfld (nx,3), STAT = istatus) ALLOCATE(rdxfld(nx,3), STAT = istatus) ALLOCATE(dyfld (ny,3), STAT = istatus) ALLOCATE(rdyfld(ny,3), STAT = istatus) CALL prepinterp(nx,ny,nxlg,nylg,nxlg_wrf,nylg_wrf,dx,dy,x,y,xs,ys,& mapproj,scalef,latnot,trulon,ctrlat,ctrlon,swx,swy, & lat_wrf,lon_wrf,x2d,y2d,iloc,jloc, & dxfld,rdxfld,dyfld,rdyfld,istatus) END IF !----------------------------------------------------------------------- ! ! Read in surface characteristic data to intialize ! !----------------------------------------------------------------------- CALL get_sfcdt(sfcinitopt(domid),sfcdtfl,ncompressx,ncompressy,nxsm,nysm,& nx,ny,nstyps,dx,dy,mapproj,trulat1,trulat2,trulon,sclfct, & ctrlat,ctrlon,soiltyp,stypfrct,vegtyp,veg,tem1, & hinterp_needed,nx_wrf,ny_wrf,dx_wrf,dy_wrf, & mapproj_wrf,trulat1_wrf,trulat2_wrf,trulon_wrf,jday, & x2d(:,:,1),y2d(:,:,1),xs,ys,iloc(:,:,1),jloc(:,:,1), & hgt_wrf,soiltyp_wrf,vegtyp_wrf,vegfrct_wrf,xice,xland, & tmn_wrf,shdmin,shdmax,albbck,snoalb,istatus) IF (istatus /= 0) & CALL arpsstop('ERROR: reading surface data. Aborted.',1) !----------------------------------------------------------------------- ! ! Set up soil model vertical layers ! !----------------------------------------------------------------------- CALL init_soil_depth(sf_surface_physics,zs_wrf,dzs_wrf,nzsoil_wrf) ! !----------------------------------------------------------------------- ! ! Get Map scale factor ! !----------------------------------------------------------------------- ! CALL lattomf(nx_wrf,ny_wrf,lat_wrf(:,:,1),msft_wrf) !mass points CALL lattomf(nx_wrf,ny_wrf,lat_wrf(:,:,2),msfu_wrf) !U points CALL lattomf(nx_wrf,ny_wrf,lat_wrf(:,:,3),msfv_wrf) !V points END IF ! ifile == 1 !======================================================================= ! ! All files should do the followings ! !======================================================================= ! Restore to WRF map projection CALL setmapr(mapproj_wrf,sclfct_wrf,lattru_wrf,lontru_wrf) CALL setorig( 1, swx_wrf, swy_wrf) !----------------------------------------------------------------------- ! ! Begin ARPS data conversions ! ! NOTE: ! From now on, pprt, ptprt, qvprt, uprt,vprt, wprt will hold total ! values of the variables instead of only the perturbation part. ! !----------------------------------------------------------------------- ! DO k=1,nz DO j=1,ny DO i=1,nx pprt(i,j,k) = pprt(i,j,k) + pbar(i,j,k) ptprt(i,j,k)= ptprt(i,j,k)+ ptbar(i,j,k) qvprt(i,j,k)= qvprt(i,j,k)+ qvbar(i,j,k) uprt(i,j,k) = uprt(i,j,k) + ubar(i,j,k) vprt(i,j,k) = vprt(i,j,k) + vbar(i,j,k) wprt(i,j,k) = wprt(i,j,k) + wbar(i,j,k) qi(i,j,k) = qi(i,j,k) + qs(i,j,k) + qh(i,j,k) END DO END DO END DO ! ! Put temperature into temp ! DO k=1,nz DO j=1,ny DO i=1,nx temp(i,j,k)=ptprt(i,j,k)*(pprt(i,j,k)/p0)**rddcp END DO END DO END DO !---------------------------------------------------------------------- ! ! Check ARPS vertical domain ! !---------------------------------------------------------------------- DO j = 1,ny-1 DO i = 1,nx-1 IF(pprt(i,j,nz-1) > ptop(domid) + 1000) THEN WRITE(6,'(/2a)') 'It seems that WRF vertical domain is ', & 'outside of the ARPS physical domain.' WRITE(6,'(/a,F6.0,a,/a,I3,a,I3,a,F6.0,a,a)') & 'WRF top pressure is ', ptop(domid),' Pascal.', & 'ARPS top pressure at i = ',i,' j = ',j, & ' is ',pprt(i,j,nz-1), ' Pascal.' CALL arpsstop('Vertical domain unmatch.',1) END IF END DO END DO ! !----------------------------------------------------------------------- ! ! Horizontally interpolate ARPS state variables to WRF grid ! !----------------------------------------------------------------------- ! IF (hinterp_needed) THEN CALL hinterp(nx,ny,nz,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & temp,dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & t_tmp,tem1,tem2,tem3) CALL hinterp(nx,ny,nz,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & pprt,dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & p_tmp,tem1,tem2,tem3) CALL hinterp(nx,ny,nz,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & zps,dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & zps_tmp,tem1,tem2,tem3) CALL hinterp(nx,ny,nz,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & qvprt,dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & qv_tmp,tem1,tem2,tem3) ELSE DO k = 1,nz DO j = 1,ny_wrf DO i = 1,nx_wrf t_tmp (i,j,k) = temp (i+1,j+1,k) p_tmp (i,j,k) = pprt (i+1,j+1,k) zps_tmp(i,j,k) = zps (i+1,j+1,k) qv_tmp (i,j,k) = qvprt(i+1,j+1,k) END DO END DO END DO END IF IF (wrftrnopt(domid) == 0 .OR. sfcinitopt(domid)(1:4) == 'ARPS') THEN IF (hinterp_needed) THEN CALL hinterp(nx,ny,1,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & zp(:,:,2),dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & hterain_wrf,tem1,tem2,tem3) ELSE DO j = 1,ny_wrf DO i = 1,nx_wrf hterain_wrf(i,j) = zp(i+1,j+1,2) ! terrain height of source data END DO END DO END IF WHERE(hterain_wrf < 0) hterain_wrf = 0.0 ELSE IF (wrftrnopt(domid) == 1) THEN hterain_wrf(:,:) = hgt_wrf(:,:) ELSE WRITE(6,'(a,I2,a)') ' **** WARNING: wrong wrftrnopt = ',wrftrnopt(domid),' ****' CALL arpsstop('Wrong wrftrnopt setting.',1) END IF ! !----------------------------------------------------------------------- ! ! Convert specific humidity to mixing ratio (kg/kg) ! ! NOTE: qv_tmp will be mixing ratio instead of specific humidity ! from now on. ! !----------------------------------------------------------------------- ! DO k = 1,nz DO j = 1,ny_wrf DO i = 1,nx_wrf qv_tmp(i,j,k) = qv_tmp(i,j,k) / (1-qv_tmp(i,j,k)) END DO END DO END DO ! !----------------------------------------------------------------------- ! ! Compute ETAP in ARPS vertical grid ! !----------------------------------------------------------------------- ! CALL compute_eta_3d(nx_wrf,ny_wrf,nz,p_tmp,t_tmp,qv_tmp,zps_tmp, & hterain_wrf,ptop(domid),etap(:,:,:,1),mu) ! T points DO k = 1, nz DO j = 1, ny_wrf DO i = 2, nx_wrf etap(i,j,k,2) = 0.5*(etap(i-1,j,k,1)+etap(i,j,k,1)) END DO etap(1,j,k,2) = 2*etap(2,j,k,2) - etap(3,j,k,2) ! U points END DO END DO DO k = 1, nz DO j = 2, ny_wrf DO i = 1, nx_wrf etap(i,j,k,3) = 0.5*(etap(i,j-1,k,1)+etap(i,j,k,1)) END DO END DO etap(:,1,k,3) = 2*etap(:,2,k,3) - etap(:,3,k,3) ! V points END DO !----------------------------------------------------------------------- ! ! Process U wind ! !----------------------------------------------------------------------- IF (hinterp_needed) THEN CALL hinterp(nx,ny,nz,nx_wrf,ny_wrf,iorder, & iloc(:,:,2),jloc(:,:,2),x,ys,x2d(:,:,2),y2d(:,:,2), & uprt,dxfld(:,2),dyfld(:,2),rdxfld(:,2),rdyfld(:,2), & work3d,tem1,tem2,tem3) ELSE DO k = 1,nz DO j = 1,ny_wrf DO i = 1,nx_wrf work3d(i,j,k) = uprt(i+1,j+1,k) END DO END DO END DO END IF CALL vinterp(nx_wrf,ny_wrf,nz,nz_wrf-1,korder,etap(:,:,:,2), & zlevels_half, work3d, u_wrf, .TRUE.) !----------------------------------------------------------------------- ! ! Process V wind ! !----------------------------------------------------------------------- IF (hinterp_needed) THEN CALL hinterp(nx,ny,nz,nx_wrf,ny_wrf,iorder, & iloc(:,:,3),jloc(:,:,3),xs,y,x2d(:,:,3),y2d(:,:,3), & vprt,dxfld(:,3),dyfld(:,3),rdxfld(:,3),rdyfld(:,3), & work3d,tem1,tem2,tem3) ELSE DO k = 1,nz DO j = 1,ny_wrf DO i = 1,nx_wrf work3d(i,j,k) = vprt(i+1,j+1,k) END DO END DO END DO END IF CALL vinterp(nx_wrf,ny_wrf,nz,nz_wrf-1,korder,etap(:,:,:,3), & zlevels_half,work3d, v_wrf, .TRUE.) !------------------------------------------------------------------- ! ! Rotate U/V from ARPS grid to WRF grid ! !------------------------------------------------------------------ IF(hinterp_needed) THEN IF (ALLOCATED(uatv_wrf)) DEALLOCATE(uatv_wrf,vatu_wrf) ALLOCATE(uatv_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) ALLOCATE(vatu_wrf (nx_wrf,ny_wrf,nz_wrf), STAT = istatus) CALL rotate_UV(nx_wrf,ny_wrf,nz_wrf, & mapproj,scalef,latnot,trulon,swx,swy,mapproj_wrf, & sclfct_wrf,lattru_wrf,lontru_wrf,swx_wrf,swy_wrf, & lon_wrf(:,:,2),lon_wrf(:,:,3),u_wrf,v_wrf, & uatv_wrf,vatu_wrf, & tem1_wrf,tem2_wrf,tem3_wrf,tem4_wrf,istatus) END IF w_wrf(:,:,:) = 0.0 !----------------------------------------------------------------------- ! ! Process Potential temperature ! !----------------------------------------------------------------------- IF (hinterp_needed) THEN CALL hinterp(nx,ny,nz,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & ptprt,dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & work3d,tem1,tem2,tem3) ELSE DO k = 1,nz DO j = 1,ny_wrf DO i = 1,nx_wrf work3d(i,j,k) = ptprt(i+1,j+1,k) END DO END DO END DO END IF CALL vinterp(nx_wrf,ny_wrf,nz,nz_wrf-1,korder,etap(:,:,:,1), & zlevels_half, work3d, pot_wrf,.TRUE.) !----------------------------------------------------------------------- ! ! Process Water vapor mixing ratio ! !----------------------------------------------------------------------- CALL vinterp(nx_wrf,ny_wrf,nz,nz_wrf-1,korder,etap(:,:,:,1), & zlevels_half,qv_tmp,qv_wrf,.TRUE.) !----------------------------------------------------------------------- ! ! Compute geopotential, air mass etc. ! ! NOTE: ! pot_wrf is potential temperature at WRF grid, total ! qv_wrf is WRF water vapor mixing ratio ! pt_wrf is the output of the purturbation potential temperature ! pot_wrf - t0 ! !----------------------------------------------------------------------- CALL compute_ph(nx_wrf,ny_wrf,nz_wrf,hterain_wrf,zlevels_wrf,ptop(domid), & mu,pot_wrf,qv_wrf,mup_wrf,mub_wrf,ph_wrf,phb_wrf, & pt_wrf,p_wrf,pb_wrf, & pt_init_wrf,tem1_wrf,tem2_wrf,tem3_wrf) !======================================================================= ! ! Input file only block ! !======================================================================= IF (ifile == 1) THEN ! compute WRF input file variables only !----------------------------------------------------------------------- ! ! Microphysics variables ! !----------------------------------------------------------------------- ! Compute QCLOUD and QRAIN, QSNOW, QICE, QGRAUP ! ! NOTE: They are not all zeros according to real.exe. ! We can initialize them here use those values from ARPS. ! IF (P_QC > 0) THEN ! QCLOUD IF (hinterp_needed) THEN CALL hinterp(nx,ny,nz,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & qc,dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & work3d,tem1,tem2,tem3) ELSE DO k = 1,nz DO j = 1,ny_wrf DO i = 1,nx_wrf work3d(i,j,k) = qc(i+1,j+1,k) END DO END DO END DO END IF CALL vinterp(nx_wrf,ny_wrf,nz,nz_wrf-1,korder,etap(:,:,:,1), & zlevels_half, work3d, qc_wrf,.TRUE.) END IF IF (P_QR > 0) THEN ! QRAIN IF (hinterp_needed) THEN CALL hinterp(nx,ny,nz,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & qr,dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & work3d,tem1,tem2,tem3) ELSE DO k = 1,nz DO j = 1,ny_wrf DO i = 1,nx_wrf work3d(i,j,k) = qr(i+1,j+1,k) END DO END DO END DO END IF CALL vinterp(nx_wrf,ny_wrf,nz,nz_wrf-1,korder,etap(:,:,:,1), & zlevels_half, work3d, qr_wrf, .TRUE.) END IF IF (P_QI > 0) THEN ! QICE IF (hinterp_needed) THEN CALL hinterp(nx,ny,nz,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & qi,dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & work3d,tem1,tem2,tem3) ELSE DO k = 1,nz DO j = 1,ny_wrf DO i = 1,nx_wrf work3d(i,j,k) = qi(i+1,j+1,k) END DO END DO END DO END IF CALL vinterp(nx_wrf,ny_wrf,nz,nz_wrf-1,korder,etap(:,:,:,1), & zlevels_half, work3d, qi_wrf, .TRUE.) END IF IF (P_QS > 0) THEN ! QSNOW IF (hinterp_needed) THEN CALL hinterp(nx,ny,nz,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & qs,dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & work3d,tem1,tem2,tem3) ELSE DO k = 1,nz DO j = 1,ny_wrf DO i = 1,nx_wrf work3d(i,j,k) = qs(i+1,j+1,k) END DO END DO END DO END IF CALL vinterp(nx_wrf,ny_wrf,nz,nz_wrf-1,korder,etap(:,:,:,1), & zlevels_half, work3d, qs_wrf,.TRUE.) END IF IF (P_QG > 0) THEN ! QGRAUP IF (hinterp_needed) THEN CALL hinterp(nx,ny,nz,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & qh,dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & work3d,tem1,tem2,tem3) ELSE DO k = 1,nz DO j = 1,ny_wrf DO i = 1,nx_wrf work3d(i,j,k) = qh(i+1,j+1,k) END DO END DO END DO END IF CALL vinterp(nx_wrf,ny_wrf,nz,nz_wrf-1,korder,etap(:,:,:,1), & zlevels_half, work3d, qg_wrf, .TRUE.) END IF !---------------------------------------------------------------------- ! ! Compute Soil related variables ! !----------------------------------------------------------------------- ! ! Interpolate soil variables ! IF (ALLOCATED(tsoil_tmp)) DEALLOCATE(tsoil_tmp,qsoil_tmp,zpsoil_tmp,tslb_wrf,smois_wrf) ALLOCATE(tsoil_tmp (nx_wrf,ny_wrf,nzsoil), STAT = istatus) ALLOCATE(qsoil_tmp (nx_wrf,ny_wrf,nzsoil), STAT = istatus) ALLOCATE(zpsoil_tmp(nx_wrf,ny_wrf,nzsoil), STAT = istatus) ALLOCATE(tslb_wrf (nx_wrf,ny_wrf,nzsoil_wrf), STAT = istatus) ALLOCATE(smois_wrf (nx_wrf,ny_wrf,nzsoil_wrf), STAT = istatus) tslb_wrf = 0.0 smois_wrf = 0.0 DO k = 1, nzsoil DO j = 1,ny DO i = 1,nx zpsoil(i,j,k) = zp(i,j,2) - zpsoil(i,j,k) END DO END DO END DO IF (hinterp_needed) THEN CALL hinterp(nx,ny,nzsoil,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & tsoil(:,:,:,0),dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & tsoil_tmp,tem1,tem2,tem3) !WDT - FIXME - be careful about horizontal interpolation of soiL !temperature and moisture without matching soil types! !See how ext2arps handles this. GMB CALL hinterp(nx,ny,nzsoil,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & qsoil(:,:,:,0),dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & qsoil_tmp,tem1,tem2,tem3) CALL hinterp(nx,ny,nzsoil,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & zpsoil,dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1),& zpsoil_tmp,tem1,tem2,tem3) ELSE DO k = 1,nzsoil DO j = 1,ny_wrf DO i = 1,nx_wrf ! tsoil_tmp(i,j,k) = tsoil(i+1,j+1,k,0) ! qsoil_tmp(i,j,k) = qsoil(i+1,j+1,k,0) !WDT 2004-01-10 GMB: for water, set WRF deep layer to be ARPS top layer !so soil depth interpolation is disabled IF ( vegtyp(i+1,j+1) == 14 .AND. k > 1) THEN tsoil_tmp(i,j,k) = tsoil(i+1,j+1,1,0) qsoil_tmp(i,j,k) = qsoil(i+1,j+1,1,0) ELSE tsoil_tmp(i,j,k) = tsoil(i+1,j+1,k,0) qsoil_tmp(i,j,k) = qsoil(i+1,j+1,k,0) ENDIF zpsoil_tmp(i,j,k) = zpsoil(i+1,j+1,k) END DO END DO END DO END IF CALL vinterp_soil(nx_wrf,ny_wrf,nzsoil,zpsoil_tmp,tsoil_tmp, & qsoil_tmp,nzsoil_wrf,zs_wrf,tslb_wrf,smois_wrf) WHERE(smois_wrf < 0.0) smois_wrf = 0.0 WHERE(smois_wrf > 1.0) smois_wrf = 1.0 IF (ANY(tslb_wrf < 170.0 .OR. tslb_wrf > 400) ) THEN WRITE(6,'(2(1x,a,/),6(10x,a,/),1x,a)') & '*********************** IMPORTANT **************************', & 'WARNING: Soil temperature is not valid.', & 'One possible reason is that flag "sfcout" is not turned on in ',& 'the ARPS data file. Please check the namelist file for flag ', & '"sfcout" and regenerate the ARPS data file. The program can ', & 'use average annual temperature read in from static file ', & 'temporarily for a replacement. But you should modify the ', & 'source code and use it at your own risk.', & '*********************** IMPORTANT **************************' CALL arpsstop('tsoil no available.',1) ! DO k = 1,nz_soil_wrf ! DO j = 1,ny_wrf ! DO i = 1,nx_wrf ! tslb_wrf(i,j,k) = tmn_wrf(i,j) ! END DO ! END DO ! END DO END IF ! ! Process Water equivalent of accumulated snow ! IF (hinterp_needed) THEN CALL hinterp(nx,ny,1,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & snowdpth,dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & snowh_wrf,tem1,tem2,tem3) ELSE snowh_wrf = snowdpth(2:nx-1,2:ny-1) END IF WHERE(snowh_wrf < 0.0) snowh_wrf = 0.0 ! ! Canopy water amount (meter, in ARPS, kg/m**2 in WRF) ! IF (hinterp_needed) THEN CALL hinterp(nx,ny,1,nx_wrf,ny_wrf,iorder, & iloc(:,:,1),jloc(:,:,1),xs,ys,x2d(:,:,1),y2d(:,:,1), & wetcanp(:,:,0), & dxfld(:,1),dyfld(:,1),rdxfld(:,1),rdyfld(:,1), & canwat_wrf,tem1,tem2,tem3) ELSE canwat_wrf(:,:) = wetcanp(2:nx-1,2:ny-1,0) END IF canwat_wrf(:,:) = canwat_wrf(:,:)*1000. ! ! SST, ARPS does not provide, use tsoil at surface to approximate ! sst_wrf(:,:) = 0.0 DO j = 1, ny_wrf-1 DO i = 1, nx_wrf-1 IF(vegtyp_wrf(i,j) == ISWATER) THEN sst_wrf(i,j) = tsoil_tmp(i,j,1) !WDT 2004-01-20 GMB: set all SST's ! Use average of SST of neighbors +-2 points away ! ! NOTE: Since WRF arrays do not provide overlay region for MPI-mode ! in current code. User should expect inconsistent with ! no-mpi mode. However, the difference should be small. ! ! No meaning to set SST over land, comment out by Y. Wang. on 02/10/2005. ! ! ELSE ! cnt = 0.125 ! sst_sum = 0.125*tsoil_tmp(i,j,nzsoil) ! own deep soil temp counts too ! DO jj = j-2,j+2 ! DO ii = i-2,i+2 ! IF (jj>0 .AND. jj< ny_wrf-1 .AND. ii>0 .AND. ii<nx_wrf-1 & ! .AND. vegtyp_wrf(ii,jj)==ISWATER) THEN ! dist_weight = 1./((jj-j)**2+(ii-i)**2) ! cnt = cnt + dist_weight ! sst_sum = sst_sum + dist_weight*tsoil_tmp(i,j,1) ! ENDIF ! END DO ! END DO ! IF (cnt > 0) sst_wrf(i,j) = sst_sum/cnt ENDIF END DO END DO !----------------------------------------------------------------------- ! ! Consistence check ! !----------------------------------------------------------------------- ! oops1 = 0 ! oops2 = 0 ! DO j = 1,ny_wrf ! DO i = 1,nx_wrf ! IF ( ( (xland(i,j) > 1.5) .AND. (vegtyp_wrf(i,j) /= ISWATER .OR. soiltyp_wrf(i,j) /= 14) ) & ! .OR. ( (xland(i,j) < 1.5) .AND. (vegtyp_wrf(i,j) == ISWATER .OR. soiltyp_wrf(i,j) == 14) ) & ! ) THEN ! IF ( sst(i,j) < 1. ) THEN ! oops1=oops1+1 ! vegtyp_wrf(i,j) = 5 ! soiltyp_wrf(i,j) = 8 ! xland(i,j) = 1 ! ELSE IF ( sst_wrf(i,j) > 1. ) THEN ! oops2=oops2+1 ! vegtyp_wrf(i,j) = ISWATER ! soiltyp_wrf(i,j) = 14 ! xland(i,j) = 2 ! ELSE ! print *,'the landmask and soil/veg cats do not match' ! print *,'i,j=',i,j ! print *,'xland=',xland(i,j) ! print *,'ivgtyp=',vegtyp_wrf(i,j) ! print *,'isltyp=',soiltyp_wrf(i,j) ! print *,'iswater=', ISWATER ! print *,'tslb=',tslb_wrf(i,j,:) ! print *,'sst=',sst_wrf(i,j) ! END IF ! END IF ! END DO ! END DO !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! OUTPUT of WRF input begin ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! !----------------------------------------------------------------------- ! ! Open WRF input file ! !----------------------------------------------------------------------- WRITE(filename,'(a,I2.2)') 'wrfinput_d',domid lenstr = LEN_TRIM(dirname) IF(lenstr > 0) filename = dirname(1:lenstr) // TRIM(filename) IF(myproc == 0) PRINT *, ' Output file name is ',TRIM(filename) CALL open_output_file(filename,'INPUT',io_form,nxlg_wrf,nylg_wrf, & nz_wrf, nzsoil_wrf,spec_bdy_width,ncid) !----------------------------------------------------------------------- ! ! Initialize and write global attributes ! !----------------------------------------------------------------------- WRITE(times_str,'(I4.4,a1,I2.2,a1,I2.2,a1,I2.2,a1,I2.2,a1,I2.2,a)')& year,'-',month,'-',day,'_',hour,':',minute,':',second,'.0000' ! set damp_opt = 0 CALL set_global_meta(nxlg_wrf,nylg_wrf,nz_wrf,execname,times_str, & domid,parent_id(domid),i_parent_start(domid), & j_parent_start(domid),grid_ratio, & dx_wrf,dy_wrf,dt,dyn_opt,diff_opt,km_opt,0, & khdif(domid),kvdif(domid),mp_physics(domid), & ra_lw_physics(domid), ra_sw_physics(domid), & sf_sfclay_physics(domid),sf_surface_physics(domid),& bl_pbl_physics(domid),cu_physics(domid), & ctrlat_wrf,ctrlon_wrf,trulat1_wrf,trulat2_wrf, & trulon_wrf, & year,jday,hour,minute,second,mapproj_wrf, & global_meta) CALL write_global_attribute(ncid,io_form,global_meta,.FALSE.) IF (io_form == IO_NET) THEN ! this call because we do not use WRF ! external IO package for NetCDF format CALL write_times_str(ncid,io_form,'Times', & times_str(1:19),times_str(1:19),ifile-1) END IF !----------------------------------------------------------------------- ! ! Data writing ! !----------------------------------------------------------------------- IF (qx_zero_out == 1) THEN qr_wrf(:,:,:) = 0.0 qi_wrf(:,:,:) = 0.0 qs_wrf(:,:,:) = 0.0 qg_wrf(:,:,:) = 0.0 END IF CALL write_wrf_input(ncid,io_form,times_str(1:19), & nx_wrf,ny_wrf,nz_wrf,nxlg_wrf,nylg_wrf, & nzsoil_wrf,spec_bdy_width,fzone_wrf,dx_wrf,dy_wrf, & zlevels_wrf, zlevels_half,ptop(domid), & mapproj_wrf,trulat1_wrf,trulat2_wrf,trulon_wrf, & ctrlat_wrf,ctrlon_wrf,lat_wrf(:,:,1),lon_wrf(:,:,1),& lat_wrf(:,:,2),lon_wrf(:,:,2),lat_wrf(:,:,3),lon_wrf(:,:,3), & lat_ll,lat_ul,lat_ur,lat_lr, & lon_ll,lon_ul,lon_ur,lon_lr, & msft_wrf,msfu_wrf,msfv_wrf,zs_wrf,dzs_wrf, & u_wrf,v_wrf,w_wrf,ph_wrf,phb_wrf,pt_wrf, & pt_init_wrf,p_wrf,pb_wrf,mup_wrf,mub_wrf,mu, & qv_wrf,qc_wrf,qr_wrf,qi_wrf,qs_wrf,qg_wrf, & soiltyp_wrf,vegtyp_wrf,vegfrct_wrf,xland, & xice,tmn_wrf,shdmax,shdmin,snoalb, & snowh_wrf,canwat_wrf,albbck,sst_wrf, & hterain_wrf,tsoil_tmp(:,:,1),tslb_wrf,smois_wrf, & tem1_wrf,tem2_wrf,tem3_wrf,tem4_wrf,istatus) !----------------------------------------------------------------------- ! ! Close WRF input file ! !----------------------------------------------------------------------- CALL close_output_file(ncid,io_form) ! close input file !write(0,*) 'input file closed' ! IF(ALLOCATED(tsoil_tmp)) DEALLOCATE(tsoil_tmp,qsoil_tmp,zpsoil_tmp) ! IF(ALLOCATED(vegtyp_wrf))DEALLOCATE(soiltyp_wrf,vegtyp_wrf,vegfrct_wrf) ! IF(ALLOCATED(xland)) DEALLOCATE(xland,xice) END IF ! ifile == 1 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ! ! WRF boundary dumping begin, block for all files ! !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% IF(create_bdy >= 1) THEN IF (spec_bdy_width > MIN(nx_wrf,ny_wrf)) THEN WRITE(6,'(a,/,a,I3,a,2(I4,a))') 'WARNING: ' // & 'Lateral boundary zone is larger than the local domain.', & ' Spec_bdy_width = ',spec_bdy_width, & ' Local domain size = ',nx_wrf,'X',ny_wrf,'.' WRITE(6,'(a,/)') 'WRF boundary file cannot be generated.' CALL arpsstop('Boundary zone too large.',1) END IF IF (ifile == 1) THEN CALL couple(nx_wrf,ny_wrf,nz_wrf,'U',mup_wrf,mub_wrf,u_wrf, & msfu_wrf,ubdy3dtemp1,work2d) CALL couple(nx_wrf,ny_wrf,nz_wrf,'V',mup_wrf,mub_wrf,v_wrf, & msfv_wrf,vbdy3dtemp1,work2d) CALL couple(nx_wrf,ny_wrf,nz_wrf,'T',mup_wrf,mub_wrf,pt_wrf, & msft_wrf,tbdy3dtemp1,work2d) CALL couple(nx_wrf,ny_wrf,nz_wrf,'H',mup_wrf,mub_wrf,ph_wrf, & msft_wrf,pbdy3dtemp1,work2d) CALL couple(nx_wrf,ny_wrf,nz_wrf,'T',mup_wrf,mub_wrf,qv_wrf, & msft_wrf,qbdy3dtemp1,work2d) mbdy2dtemp1 = mup_wrf ELSE CALL couple(nx_wrf,ny_wrf,nz_wrf,'U',mup_wrf,mub_wrf,u_wrf, & msfu_wrf,ubdy3dtemp2,work2d) CALL couple(nx_wrf,ny_wrf,nz_wrf,'V',mup_wrf,mub_wrf,v_wrf, & msfv_wrf,vbdy3dtemp2,work2d) CALL couple(nx_wrf,ny_wrf,nz_wrf,'T',mup_wrf,mub_wrf,pt_wrf, & msft_wrf,tbdy3dtemp2,work2d) CALL couple(nx_wrf,ny_wrf,nz_wrf,'H',mup_wrf,mub_wrf,ph_wrf, & msft_wrf,pbdy3dtemp2,work2d) CALL couple(nx_wrf,ny_wrf,nz_wrf,'T',mup_wrf,mub_wrf,qv_wrf, & msft_wrf,qbdy3dtemp2,work2d) mbdy2dtemp2 = mup_wrf count_in_loop = 0 DO WHILE (abstimec < abstimen) count_bdy = count_bdy + 1 ! count in the output boundary file count_in_loop = count_in_loop + 1 ! count in this WHILE loop abstimecn = abstimec + tintv_bdywrf CALL abss2ctim(abstimecn, year, month, day, hour, minute, second ) WRITE(nextbdytimes,'(I4.4,a1,I2.2,a1,I2.2,a1,I2.2,a1,I2.2,a1,I2.2,a)')& year,'-',month,'-',day,'_',hour,':',minute,':',second,'.0000' IF (abstimecn <= abstimep .OR. abstimecn > abstimen ) THEN WRITE(0,'(1x,a,I4,/,8x,2a,/,8x,a,I12,/,8x,a,I12,a)') & 'ERROR: data time inconsistent for boundary ',count_bdy, & 'Expecting data time at or beyond: ',nextbdytimes, & 'with abstime = ',abstimecn, & 'but found data at time = ',abstime,'.' CALL arpsstop('Data time inconsistent.',1) END IF CALL abss2ctim(abstimec, year, month, day, hour, minute, second ) CALL julday( year, month, day, jday) WRITE(times_str,'(I4.4,a1,I2.2,a1,I2.2,a1,I2.2,a1,I2.2,a1,I2.2,a)')& year,'-',month,'-',day,'_',hour,':',minute,':',second,'.0000' IF ( count_bdy == 1) THEN ! IF( ABS(abstime - abstime1 - tintv_bdyin) > a_small_real_number) THEN ! WRITE(6,'(2a)') 'ERROR: the first boundary file must be ', & ! 'data valid at: the initial time + tintv_bdyin.' ! WRITE(6,'(a,I4.4,5(a,I2.2),a,I6)') 'Expected data time: ', & ! year1,'/',month1,'/',day1,'-',hour1,':',minute1,':', & ! second1,' + ',tintv_bdyin ! WRITE(6,'(a,I4.4,5(a,I2.2))') 'Found data valid time: ', & ! year,'/',month,'/',day,'-',hour,':',minute,':',second ! CALL arpsstop('Inconsistent time.',1) ! END IF filename = 'wrfbdy_d01' lenstr = LEN_TRIM(dirname) IF(lenstr > 0) filename = dirname(1:lenstr) // TRIM(filename) IF (myproc == 0) PRINT *, ' Output file name is ',TRIM(filename) CALL open_output_file(filename,'BDY',io_form, & nxlg_wrf,nylg_wrf,nz_wrf,nzsoil_wrf,spec_bdy_width, & ncid) ! ! Initialize and write global attributes ! CALL set_global_meta(nxlg_wrf,nylg_wrf,nz_wrf,execname, & times_str,1,1,1,1,1, & dx_wrf,dy_wrf,dt,dyn_opt,diff_opt,km_opt,0, & khdif(1),kvdif(1), & mp_physics(1),ra_lw_physics(1), ra_sw_physics(1),& sf_sfclay_physics(1),sf_surface_physics(1), & bl_pbl_physics(1),cu_physics(1), & ctrlat_wrf,ctrlon_wrf,trulat1_wrf,trulat2_wrf, & trulon_wrf, & year,jday,hour,minute,second,mapproj_wrf, & global_meta) ! set damp_opt = 0 CALL write_global_attribute(ncid,io_form,global_meta,.TRUE.) END IF ! initilize boundary file IF (myproc == 0) WRITE(6,'(a,I4.4,5(a,I2.2))') & '=== Begin boundary dumps at: ', year,'/',month,'/',day, & '-',hour,':',minute,':',second IF (io_form == IO_NET) THEN ! this call because we do not use WRF ! external IO package for NetCDF format CALL write_times_str(ncid,io_form,'Times', & times_str(1:19),times_str(1:19),count_bdy) END IF CALL write_times_str(ncid,io_form,'THISBDYTIME', & times_str(1:19),times_str(1:19),count_bdy) CALL write_times_str(ncid,io_form,'NEXTBDYTIME', & times_str(1:19),nextbdytimes(1:19),count_bdy) !----------------------------------------------------------------------- ! ! Prepare boundary data in two time levels, abstimec & abstimecn ! !----------------------------------------------------------------------- IF (count_in_loop == 1) THEN !write(0,*) 'Count ',count_bdy,count_in_loop,' tempc => temp1' ubdy3dtempc => ubdy3dtemp1 vbdy3dtempc => vbdy3dtemp1 tbdy3dtempc => tbdy3dtemp1 pbdy3dtempc => pbdy3dtemp1 qbdy3dtempc => qbdy3dtemp1 mbdy2dtempc => mbdy2dtemp1 ELSE !write(0,*) 'Count ',count_bdy,count_in_loop,' tempc => tempcn' ubdy3dtempc => ubdy3dtempcn vbdy3dtempc => vbdy3dtempcn tbdy3dtempc => tbdy3dtempcn pbdy3dtempc => pbdy3dtempcn qbdy3dtempc => qbdy3dtempcn mbdy2dtempc => mbdy2dtempcn END IF IF (abstimecn == abstimen) THEN !write(0,*) 'Count ',count_bdy,count_in_loop,' tempcn => temp2' ubdy3dtempcn => ubdy3dtemp2 vbdy3dtempcn => vbdy3dtemp2 tbdy3dtempcn => tbdy3dtemp2 pbdy3dtempcn => pbdy3dtemp2 qbdy3dtempcn => qbdy3dtemp2 mbdy2dtempcn => mbdy2dtemp2 ELSE ! time interpolation - linear abstdiff1 = abstimecn - abstimep abstdiff2 = abstimen - abstimecn abstdiff = abstimen - abstimep IF (abstdiff / tintv_bdyin > 2) THEN WRITE(6,'(/,1x,a,I8,a,/,9x,a)') & 'WARNING: The time gap is too large (',abstdiff,' seconds)',& 'Time interpolation may be not accurate enough.' END IF !----------------------------------------------------------------------- ! ! This block is to make sure that we will allocate at most 2 sets of ! temporary arrays, bdy3d?temp1-5 and bdy2d?temp1 (where ? is 1/2). ! ! case 1. Only do 1 loop, no temporary arrays should be allocated (previous behavior); ! case 2. Do 2 loops, 1 set of temporary arrays should be allocated; ! case 3. Do 3 or more loops, 2 sets of temporary arrays will be allocated. ! !----------------------------------------------------------------------- IF (mod(count_in_loop,2) == 0) THEN IF (.NOT. ALLOCATED(bdy3d2temp1)) & ALLOCATE(bdy3d2temp1(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) IF (.NOT. ALLOCATED(bdy3d2temp2)) & ALLOCATE(bdy3d2temp2(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) IF (.NOT. ALLOCATED(bdy3d2temp3)) & ALLOCATE(bdy3d2temp3(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) IF (.NOT. ALLOCATED(bdy3d2temp4)) & ALLOCATE(bdy3d2temp4(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) IF (.NOT. ALLOCATED(bdy3d2temp5)) & ALLOCATE(bdy3d2temp5(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) IF (.NOT. ALLOCATED(bdy2d2temp1)) & ALLOCATE(bdy2d2temp1(nx_wrf,ny_wrf), STAT = istatus) !write(0,*) 'Count ',count_bdy,count_in_loop,' tempcn => bdy3d2' ubdy3dtempcn => bdy3d2temp1 vbdy3dtempcn => bdy3d2temp2 tbdy3dtempcn => bdy3d2temp3 pbdy3dtempcn => bdy3d2temp4 qbdy3dtempcn => bdy3d2temp5 mbdy2dtempcn => bdy2d2temp1 ELSE IF (.NOT. ALLOCATED(bdy3d1temp1)) & ALLOCATE(bdy3d1temp1(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) IF (.NOT. ALLOCATED(bdy3d1temp2)) & ALLOCATE(bdy3d1temp2(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) IF (.NOT. ALLOCATED(bdy3d1temp3)) & ALLOCATE(bdy3d1temp3(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) IF (.NOT. ALLOCATED(bdy3d1temp4)) & ALLOCATE(bdy3d1temp4(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) IF (.NOT. ALLOCATED(bdy3d1temp5)) & ALLOCATE(bdy3d1temp5(nx_wrf,ny_wrf,nz_wrf), STAT = istatus) IF (.NOT. ALLOCATED(bdy2d1temp1)) & ALLOCATE(bdy2d1temp1(nx_wrf,ny_wrf), STAT = istatus) !write(0,*) 'Count ',count_bdy,count_in_loop,' tempcn => bdy3d1' ubdy3dtempcn => bdy3d1temp1 vbdy3dtempcn => bdy3d1temp2 tbdy3dtempcn => bdy3d1temp3 pbdy3dtempcn => bdy3d1temp4 qbdy3dtempcn => bdy3d1temp5 mbdy2dtempcn => bdy2d1temp1 END IF w2 = abstdiff1*1.0/abstdiff w1 = abstdiff2*1.0/abstdiff !write(0,*) 'Count ',count_bdy,count_in_loop,' w1,w2 = ',w1,w2 ubdy3dtempcn(:,:,:) = w1*ubdy3dtemp1(:,:,:) + w2*ubdy3dtemp2(:,:,:) vbdy3dtempcn(:,:,:) = w1*vbdy3dtemp1(:,:,:) + w2*vbdy3dtemp2(:,:,:) tbdy3dtempcn(:,:,:) = w1*tbdy3dtemp1(:,:,:) + w2*tbdy3dtemp2(:,:,:) pbdy3dtempcn(:,:,:) = w1*pbdy3dtemp1(:,:,:) + w2*pbdy3dtemp2(:,:,:) qbdy3dtempcn(:,:,:) = w1*qbdy3dtemp1(:,:,:) + w2*qbdy3dtemp2(:,:,:) mbdy2dtempcn(:,:) = w1*mbdy2dtemp1(:,:) + w2*mbdy2dtemp2(:,:) END IF !----------------------------------------------------------------------- ! ! Boundary data writing ! !----------------------------------------------------------------------- CALL write_wrf_bdy(ncid,io_form,times_str(1:19),count_bdy, & nx_wrf,ny_wrf,nz_wrf,spec_bdy_width, & nxlg_wrf,nylg_wrf,fzone_wrf,tintv_bdywrf, & ubdy3dtempc,ubdy3dtempcn, & vbdy3dtempc,vbdy3dtempcn, & tbdy3dtempc,tbdy3dtempcn, & pbdy3dtempc,pbdy3dtempcn, & qbdy3dtempc,qbdy3dtempcn, & mbdy2dtempc,mbdy2dtempcn, & bdys,bdyn,bdyw,bdye, & blgw,blge,blgs,blgn, & tem1_wrf,tem2_wrf,istatus) abstimec = abstimecn END DO ! while loop ubdy3dtemp1 = ubdy3dtemp2 vbdy3dtemp1 = vbdy3dtemp2 tbdy3dtemp1 = tbdy3dtemp2 pbdy3dtemp1 = pbdy3dtemp2 qbdy3dtemp1 = qbdy3dtemp2 mbdy2dtemp1 = mbdy2dtemp2 ! close boundary file, The last ARPS file should not be skipped IF(ifile == nhisfile) THEN CALL close_output_file(ncid,io_form) !write(0,*) 'deallocating bdy arrays' ! IF (ALLOCATED(ubdy3dtemp1)) THEN ! DEALLOCATE(ubdy3dtemp1,vbdy3dtemp1,tbdy3dtemp1) ! DEALLOCATE(pbdy3dtemp1,qbdy3dtemp1,mbdy2dtemp1) ! DEALLOCATE(ubdy3dtemp2,vbdy3dtemp2,tbdy3dtemp2) ! DEALLOCATE(pbdy3dtemp2,qbdy3dtemp1,mbdy2dtemp2) ! DEALLOCATE(bdys,bdyn,bdyw,bdye) ! DEALLOCATE(blgs,blgn,blgw,blge) ! END IF END IF END IF ! ifile > 1 END IF ! create_wrf_bdy IF(ifile == nhisfile .AND. domid == 1) THEN year2 = year month2 = month day2 = day hour2 = hour minute2= minute second2= second END IF END DO ! ifile CALL io_shutdown(io_form) !----------------------------------------------------------------------- ! ! Generate a ready file if needed ! !----------------------------------------------------------------------- IF ( readyfl == 1 .AND. myproc == 0) THEN CALL getunit( nchout ) IF(create_bdy >= 1) THEN WRITE(filename,'(a,I2.2,3a)') 'wrfinput_bdy_d',domid,'.',fmtstr(io_form),'_ready' ELSE WRITE (filename,'(a,I2.2,3a)') 'wrfinput_d',domid,'.',fmtstr(io_form),'_ready' END IF WRITE(filename,'(a)') TRIM(dirname) // TRIM(filename) OPEN (UNIT=nchout,FILE=trim(filename)) WRITE (nchout,'(a,I2.2)') 'wrfinput_d',domid IF(create_bdy >= 1) WRITE(nchout,'(a)') 'wrfbdy_d01' CLOSE (nchout) CALL retunit (nchout ) END IF ELSE ! Prepare for namelist.input only nx_wrf = nx_wrfnm(domid) ny_wrf = ny_wrfnm(domid) dx_wrf = dx_wrfnm(domid) dy_wrf = dy_wrfnm(domid) CALL get_check_grid_ratio(domid,nx_wrf,ny_wrf,dx_wrf,dy_wrf, & dx_wrfnm(parent_id(domid)),dy_wrfnm(parent_id(domid)), & a_small_real_number,grid_ratio,istatus) IF (istatus /= 0) CALL arpsstop('Unacceptable parent_grid_ratio or nesting grid size.',1) parent_grid_ratio(domid) = grid_ratio year1(domid) = year1(parent_id(domid)) month1(domid) = month1(parent_id(domid)) day1(domid) = day1(parent_id(domid)) hour1(domid) = hour1(parent_id(domid)) minute1(domid) = minute1(parent_id(domid)) second1(domid) = second1(parent_id(domid)) END IF ! input_from_file END DO ! domid !----------------------------------------------------------------- ! ! Write out a namelist file for WRF to run if needed. ! !----------------------------------------------------------------- IF (create_namelist == 1 .AND. myproc == 0) THEN WRITE(6,'(/,1x,a)') 'Writing namelist input file for the WRF model ...' CALL write_wrf_namelist(nchout,TRIM(dirname)//TRIM(wrfnamelist),io_form, & max_dom,input_from_file,nx_wrfnm,ny_wrfnm,nz_wrf,nzsoil_wrf,& nproc_x,nproc_y,fzone_wrf,dx_wrfnm,dy_wrfnm, & dt,spec_bdy_width,parent_time_step_ratio, & i_parent_start,j_parent_start,parent_id,parent_grid_ratio,& year2,month2,day2,hour2,minute2,second2, tintv_bdywrf, & year1,month1,day1,hour1,minute1,second1, & mp_physics,ra_lw_physics,ra_sw_physics,sf_sfclay_physics, & sf_surface_physics,bl_pbl_physics,cu_physics, & diff_opt,km_opt,khdif,kvdif,nprocx_wrf,nprocy_wrf, & frames_per_outfile,restart_interval,radt,cudt, & ifsnow,w_damping,io_form_history,io_form_restart, & history_interval, & indir,outdir,staticdir,readyfl,pd_moist,istatus) END IF !----------------------------------------------------------------------- ! ! Deallocate arrays, not necessary. For testing purpose only ! !----------------------------------------------------------------------- !write(0,*) 'Deallocating' ! IF (ASSOCIATED(x)) THEN ! DEALLOCATE(x,y,z,zp,zpsoil) ! DEALLOCATE(uprt,vprt,wprt,ptprt,pprt,qvprt) ! DEALLOCATE(qv,qc,qr,qi,qs,qh) ! DEALLOCATE(ubar,vbar,wbar,ptbar,pbar,rhobar,qvbar) ! DEALLOCATE(soiltyp,stypfrct,vegtyp,veg) ! DEALLOCATE(tsoil,qsoil,wetcanp,snowdpth) ! DEALLOCATE(dum2da, dum3da) ! END IF ! ! IF (ALLOCATED(xs)) DEALLOCATE(xs,ys,zps,temp) ! ! IF (ALLOCATED(zlevels_half)) THEN !write(0,*) 'deallocating zlevels_half' ! DEALLOCATE(zlevels_half) !write(0,*) 'deallocating mu' ! DEALLOCATE(mu,hterain_wrf) !write(0,*) 'deallocating integer' ! DEALLOCATE(soiltyp_wrf,vegtyp_wrf) !write(0,*) 'deallocating next' ! DEALLOCATE(vegfrct_wrf,xice,xland) !write(0,*) 'deallocating sst_wrf' ! DEALLOCATE(sst_wrf,hgt_wrf,tmn_wrf,shdmin,shdmax,albbck,snoalb) !write(0,*) 'deallocating snowh_wrf' ! DEALLOCATE(snowh_wrf,canwat_wrf) ! END IF ! ! IF (ALLOCATED(u_wrf)) THEN ! DEALLOCATE(u_wrf,v_wrf,w_wrf,ph_wrf,phb_wrf,pot_wrf) ! DEALLOCATE(pt_wrf,p_wrf,pb_wrf,qv_wrf,qc_wrf) ! DEALLOCATE(qr_wrf,qi_wrf,qs_wrf,qg_wrf,mup_wrf,mub_wrf) ! DEALLOCATE(tem1_wrf,tem2_wrf,tem3_wrf,tem4_wrf) ! DEALLOCATE(zs_wrf,dzs_wrf,msft_wrf,msfu_wrf,msfv_wrf) ! END IF ! ! IF (ALLOCATED(lat_wrf)) DEALLOCATE(lat_wrf,lon_wrf) ! ! IF (ALLOCATED(tem1)) THEN ! DEALLOCATE(tem1,tem2,tem3) ! DEALLOCATE(t_tmp,p_tmp,qv_tmp,zps_tmp) ! DEALLOCATE(etap,work2d,work3d) ! END IF ! ! IF (ALLOCATED(x2d)) DEALLOCATE(x2d,y2d,iloc,jloc) ! !! IF (ASSOCIATED(xsm)) THEN !! DEALLOCATE(xsm,ysm,zpsm,zpsoilsm) !! DEALLOCATE(uprtsm,vprtsm,wprtsm,ptprtsm,pprtsm,qvprtsm) !! DEALLOCATE(qcsm,qrsm,qism,qssm,qhsm) !! DEALLOCATE(ubarsm,vbarsm,wbarsm,ptbarsm,pbarsm,rhobarsm,qvbarsm) !! DEALLOCATE(soiltypsm,stypfrctsm,vegtypsm,vegsm) !! DEALLOCATE(tsoilsm,qsoilsm,wetcanpsm,snowdpthsm) !! END IF ! ! IF (ALLOCATED(dxfld)) DEALLOCATE(dxfld,dyfld,rdxfld,rdyfld) ! ! IF (ALLOCATED(uatv_wrf)) DEALLOCATE(uatv_wrf,vatu_wrf) ! ! IF (ALLOCATED(tsoil_tmp)) DEALLOCATE(tsoil_tmp,qsoil_tmp,zpsoil_tmp,tslb_wrf,smois_wrf) ! IF (mp_opt > 0) THEN IF (myproc == 0) WRITE(6,'(/,5x,a)') '==== ARPS2WRF_MPI terminated normally ====' CALL mpexit(0) ELSE WRITE(6,'(/,5x,a)') '==== ARPS2WRF terminated normally ====' STOP END IF END PROGRAM arps2wrf