PROGRAM arpsensic,110 ! !################################################################## !################################################################## !###### ###### !###### PROGRAM ARPSENSIC ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! ! !----------------------------------------------------------------------- ! ! PURPOSE: ! Generate one perturbation initial condition files from two ! sets of ARPS output files and write it out for the use in ! ENSEMBLE forecast. ! The idea is based on the SLAF (Scaled Lagged Average Forecast). ! The procedure is as follows: ! 1, read file a; ! 2, read in file b; ! 3, find the difference bwtween a and b, store it in a. a=a-b ! 4, read in the base file c or use b as the control file (iread=0). ! 5, generate the perturbation c=c+a*n (b is used as c in code) ! 6, n is input as the real variable iorder, a real factor ! ! It shares with the model the include file dims.inc for ! definition of dimensions and domain size. a,b,c have the same ! dimensions and the same grid structure. ! ! Parameters grdout,varout,mstout,iceout and trbout should be input ! with the same values as in the data dump subroutines in the model. ! ! AUTHOR: Dingchen Hou ! History: Apr. 30, 1998: developed from the framework of ARPSDIFF. ! Sep. 15, 1999: modified to include soil variables in ! perturbation change input to namelist format. ! Feb-Apr, 2002: (F.KONG) Major modifications to: ! - accommodate BGM (Breeding Fast Growing Mode) IC ! generation, including the initial random perturbation ! procedure (with inibred = 1 and specified iseed). ! During the regular breeding cycles, certain scale ! factors are calculated to control the amplitude of ! the growing perturbations (with iensopt = 1). When ! generating BGM IC, the two 12h forecasts from the ! paired breeding members are specified as a and b, ! and the analysis data (c) must be read in (iread=1). ! ! When generating initial BGM IC, the only one data ! needed is the analysis (both a and b) ! ! - read in domain config directly from data ! ! MODIFIED: ! ! 05/28/2002 (J. Brotzge) ! Added tsoil/qsoil to accomodate new soil scheme. ! ! 1 June 2002 Eric Kemp ! Soil variable updates ! ! 02/01/2005 (F. KONG) ! Change iorder to real factor (originally interger) ! ! 04/10/2005 (F. KONG) ! Change random perturbation to normal distribution ! with zero mean and specified variance (add subroutine ! normalrand for this purpose). Specify variance via ! arpsensic.input; Add computer assigning seed value ! ! 02/15/2007 (F. KONG) ! Major rewrite in this revision to: ! - Changed to MPI version. ! - Remove domain matching check to streamline the code ! (Always assuming all read in data have the same ARPS ! domain seeting) ! !----------------------------------------------------------------------- ! ! DATA ARRAYS READ IN: ! ! x x coordinate of grid points in computational space (m) ! y y coordinate of grid points in computational 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 grid points in the soil (m) ! ! uprt perturbation x component of velocity (m/s) ! vprt perturbation y component of velocity (m/s) ! wprt perturbation z component of velocity (m/s) ! ptprt perturbation potential temperature (K) ! pprt perturbation pressure (Pascal) ! qvprt perturbation 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) ! ! tke ! kmh ! kmv ! ! soiltyp Soil type ! stypfrct ! vegtyp Vegetation type ! lai Leaf Area Index ! roufns Surface roughness ! veg Vegetation fraction ! ! tsoil Soil temperature (K) ! qsoil Soil moisture (m**3/m**3) ! wetcanp Canopy water amount ! snowdpth ! ! raing Grid supersaturation rain ! rainc Cumulus convective rain ! ! ! WORK ARRAYS: ! ! tem1 Temporary work array. ! tem2 Temporary work array. ! tem3dsoil Work array ! vtem3dsoil Work array !----------------------------------------------------------------------- ! ! Variable Declarations: ! !----------------------------------------------------------------------- ! IMPLICIT NONE ! !----------------------------------------------------------------------- ! ! Include files: ! !----------------------------------------------------------------------- ! INCLUDE 'globcst.inc' INCLUDE 'grid.inc' INCLUDE 'indtflg.inc' INCLUDE 'mp.inc' INTEGER :: nx,ny,nz,nzsoil INTEGER :: nstyps ! Maximum number of soil types in each ! grid box ! !----------------------------------------------------------------------- ! ! Arrays to be read in: ! !----------------------------------------------------------------------- ! REAL, ALLOCATABLE :: x (:) ! The x-coord. of the physical and ! computational grid. Defined at u-point. REAL, ALLOCATABLE :: y (:) ! The y-coord. of the physical and ! computational grid. Defined at v-point. REAL, ALLOCATABLE :: z (:) ! The z-coord. of the computational grid. ! Defined at w-point on the staggered grid. REAL, ALLOCATABLE :: zp (:,:,:) ! The physical height coordinate defined at ! w-point of the staggered grid. REAL, ALLOCATABLE :: zpsoil(:,:,:) ! The physical height coordinate defined at ! w-point of the soil grid. REAL, ALLOCATABLE :: uprt (:,:,:) ! Perturbation u-velocity (m/s) REAL, ALLOCATABLE :: vprt (:,:,:) ! Perturbation v-velocity (m/s) REAL, ALLOCATABLE :: wprt (:,:,:) ! Perturbation w-velocity (m/s) REAL, ALLOCATABLE :: ptprt (:,:,:) ! Perturbation potential temperature (K) REAL, ALLOCATABLE :: pprt (:,:,:) ! Perturbation pressure (Pascal) REAL, ALLOCATABLE :: qvprt (:,:,:) ! Perturbation water vapor specific ! humidity (kg/kg) REAL, ALLOCATABLE :: qc (:,:,:) ! Cloud water mixing ratio (kg/kg) REAL, ALLOCATABLE :: qr (:,:,:) ! Rain water mixing ratio (kg/kg) REAL, ALLOCATABLE :: qi (:,:,:) ! Cloud ice mixing ratio (kg/kg) REAL, ALLOCATABLE :: qs (:,:,:) ! Snow mixing ratio (kg/kg) REAL, ALLOCATABLE :: qh (:,:,:) ! Hail mixing ratio (kg/kg) REAL, ALLOCATABLE :: tke (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2) REAL, ALLOCATABLE :: kmh (:,:,:) ! Horizontal turb. mixing coef. for ! momentum. ( m**2/s ) REAL, ALLOCATABLE :: kmv (:,:,:) ! Vertical turb. mixing coef. for ! momentum. ( m**2/s ) REAL, ALLOCATABLE :: ubar (:,:,:) ! Base state u-velocity (m/s) REAL, ALLOCATABLE :: vbar (:,:,:) ! Base state v-velocity (m/s) REAL, ALLOCATABLE :: wbar (:,:,:) ! Base state w-velocity (m/s) REAL, ALLOCATABLE :: ptbar (:,:,:) ! Base state potential temperature (K) REAL, ALLOCATABLE :: pbar (:,:,:) ! Base state pressure (Pascal) REAL, ALLOCATABLE :: rhobar (:,:,:) ! Base state air density (kg/m**3) REAL, ALLOCATABLE :: qvbar (:,:,:) ! Base state water vapor specific ! humidity (kg/kg) INTEGER, ALLOCATABLE :: soiltyp (:,:,:) ! Soil type REAL, ALLOCATABLE :: stypfrct(:,:,:) ! Soil fraction INTEGER, ALLOCATABLE :: vegtyp (:,:) ! Vegetation type REAL, ALLOCATABLE :: lai (:,:) ! Leaf Area Index REAL, ALLOCATABLE :: roufns (:,:) ! Surface roughness REAL, ALLOCATABLE :: veg (:,:) ! Vegetation fraction REAL, ALLOCATABLE :: tsoil (:,:,:,:) ! Soil temperature (K) REAL, ALLOCATABLE :: qsoil (:,:,:,:) ! Soil moisture (m**3/m**3) REAL, ALLOCATABLE :: wetcanp(:,:,:) ! Canopy water amount REAL, ALLOCATABLE :: snowdpth(:,:) ! Snow depth (m) REAL, ALLOCATABLE :: raing(:,:) ! Grid supersaturation rain REAL, ALLOCATABLE :: rainc(:,:) ! Cumulus convective rain REAL, ALLOCATABLE :: prcrate(:,:,:) ! precipitation rate (kg/(m**2*s)) ! prcrate(1,1,1) = total precip. rate ! prcrate(1,1,2) = grid scale precip. rate ! prcrate(1,1,3) = cumulus precip. rate ! prcrate(1,1,4) = microphysics precip. rate REAL, ALLOCATABLE :: radfrc(:,:,:) ! Radiation forcing (K/s) REAL, ALLOCATABLE :: radsw (:,:) ! Solar radiation reaching the surface REAL, ALLOCATABLE :: rnflx (:,:) ! Net radiation flux absorbed by surface REAL, ALLOCATABLE :: radswnet(:,:) ! Net shortwave radiation REAL, ALLOCATABLE :: radlwin(:,:) ! Incoming longwave radiation REAL, ALLOCATABLE :: usflx (:,:) ! Surface flux of u-momentum (kg/(m*s**2)) REAL, ALLOCATABLE :: vsflx (:,:) ! Surface flux of v-momentum (kg/(m*s**2)) REAL, ALLOCATABLE :: ptsflx(:,:) ! Surface heat flux (K*kg/(m*s**2)) REAL, ALLOCATABLE :: qvsflx(:,:) ! Surface moisture flux (kg/(m**2*s) ! !----------------------------------------------------------------------- ! ! Verification Arrays ! !----------------------------------------------------------------------- ! REAL, ALLOCATABLE :: vx (:) ! The x-coord. of the physical and ! computational grid. Defined at u-point. REAL, ALLOCATABLE :: vy (:) ! The y-coord. of the physical and ! computational grid. Defined at v-point. REAL, ALLOCATABLE :: vz (:) ! The z-coord. of the computational grid. ! Defined at w-point on the staggered grid. REAL, ALLOCATABLE :: vzp (:,:,:) ! The physical height coordinate defined at ! w-point of the staggered grid. REAL, ALLOCATABLE :: vzpsoil(:,:,:) ! The physical height coordinate defined at ! w-point of the soil grid. REAL, ALLOCATABLE :: vuprt (:,:,:) ! Perturbation u-velocity (m/s) REAL, ALLOCATABLE :: vvprt (:,:,:) ! Perturbation v-velocity (m/s) REAL, ALLOCATABLE :: vwprt (:,:,:) ! Perturbation w-velocity (m/s) REAL, ALLOCATABLE :: vptprt (:,:,:) ! Perturbation potential temperature (K) REAL, ALLOCATABLE :: vpprt (:,:,:) ! Perturbation pressure (Pascal) REAL, ALLOCATABLE :: vqvprt (:,:,:) ! Perturbation water vapor specific ! humidity (kg/kg) REAL, ALLOCATABLE :: vqc (:,:,:) ! Cloud water mixing ratio (kg/kg) REAL, ALLOCATABLE :: vqr (:,:,:) ! Rain water mixing ratio (kg/kg) REAL, ALLOCATABLE :: vqi (:,:,:) ! Cloud ice mixing ratio (kg/kg) REAL, ALLOCATABLE :: vqs (:,:,:) ! Snow mixing ratio (kg/kg) REAL, ALLOCATABLE :: vqh (:,:,:) ! Hail mixing ratio (kg/kg) REAL, ALLOCATABLE :: vtke (:,:,:) ! Turbulent Kinetic Energy ((m/s)**2) REAL, ALLOCATABLE :: vkmh (:,:,:) ! Horizontal turb. mixing coef. for ! momentum. ( m**2/s ) REAL, ALLOCATABLE :: vkmv (:,:,:) ! Vertical turb. mixing coef. for ! momentum. ( m**2/s ) REAL, ALLOCATABLE :: vubar (:,:,:) ! Base state u-velocity (m/s) REAL, ALLOCATABLE :: vvbar (:,:,:) ! Base state v-velocity (m/s) REAL, ALLOCATABLE :: vwbar (:,:,:) ! Base state w-velocity (m/s) REAL, ALLOCATABLE :: vptbar (:,:,:) ! Base state potential temperature (K) REAL, ALLOCATABLE :: vpbar (:,:,:) ! Base state pressure (Pascal) REAL, ALLOCATABLE :: vrhobar (:,:,:) ! Base state air density (kg/m**3) REAL, ALLOCATABLE :: vqvbar (:,:,:) ! Base state water vapor specific ! humidity (kg/kg) INTEGER, ALLOCATABLE :: vsoiltyp (:,:,:) ! Soil type REAL, ALLOCATABLE :: vstypfrct(:,:,:) ! Soil fraction INTEGER, ALLOCATABLE :: vvegtyp (:,:) ! Vegetation type REAL, ALLOCATABLE :: vlai (:,:) ! Leaf Area Index REAL, ALLOCATABLE :: vroufns (:,:) ! Surface roughness REAL, ALLOCATABLE :: vveg (:,:) ! Vegetation fraction REAL, ALLOCATABLE :: vtsoil (:,:,:,:) ! Soil temperature (K) REAL, ALLOCATABLE :: vqsoil (:,:,:,:) ! Soil moisture (m**3/m**3) REAL, ALLOCATABLE :: vwetcanp(:,:,:) ! Canopy water amount REAL, ALLOCATABLE :: vsnowdpth(:,:) ! Snow depth (m) REAL, ALLOCATABLE :: vraing(:,:) ! Grid supersaturation rain REAL, ALLOCATABLE :: vrainc(:,:) ! Cumulus convective rain ! !----------------------------------------------------------------------- ! ! Work Arrays ! !----------------------------------------------------------------------- ! REAL, ALLOCATABLE :: tem1(:,:,:) REAL, ALLOCATABLE :: tem2(:,:,:) REAL, ALLOCATABLE :: tem3(:,:,:) REAL, ALLOCATABLE :: tem3dsoil(:,:,:) REAL, ALLOCATABLE :: vtem1(:,:,:) REAL, ALLOCATABLE :: vtem2(:,:,:) REAL, ALLOCATABLE :: vtem3(:,:,:) REAL, ALLOCATABLE :: vtem3dsoil(:,:,:) REAL, ALLOCATABLE :: xs(:) REAL, ALLOCATABLE :: ys(:) REAL, ALLOCATABLE :: zps(:,:,:) REAL, ALLOCATABLE :: lat(:,:),lon(:,:) REAL, ALLOCATABLE :: vxs(:) REAL, ALLOCATABLE :: vys(:) REAL, ALLOCATABLE :: vzps(:,:,:) REAL, ALLOCATABLE :: dxfld(:) REAL, ALLOCATABLE :: dyfld(:) REAL, ALLOCATABLE :: rdxfld(:) REAL, ALLOCATABLE :: rdyfld(:) ! !----------------------------------------------------------------------- ! ! Misc. local variables ! !----------------------------------------------------------------------- ! CHARACTER (LEN=80) :: fcrnam,runnmin CHARACTER (LEN=256) :: filename(3),grdbasfn(3) CHARACTER (LEN=256) :: fbasfn,filnam CHARACTER (LEN=40) :: q3dname,q3dunit INTEGER :: lengbf(3),lenfil(3) INTEGER :: ifproj,ivproj REAL :: flatnot(2),vlatnot(2) REAL :: fscale,ftrulon,fdx,fdy,fx0,fy0 REAL :: fctrlat,fctrlon REAL :: vscale,vtrulon,vdx,vdy,vx0,vy0 REAL :: vctrlat,vctrlon REAL :: time,xctr,yctr INTEGER :: i,j,k INTEGER :: grdbas INTEGER :: hinfmt,iread,isread,iensopt,inibred,iseed REAL :: iorder,uvar,vvar,wvar,ptvar,qvvar INTEGER :: ireturn LOGICAL :: comcoord INTEGER :: nchin INTEGER :: zpsoilin,tsoilin,qsoilin,wcanpin,snowdin INTEGER :: dmp_out_joined, pertout ! INTEGER :: istatus INTEGER :: idate(8) REAL :: utot,utot2,usd,uscl REAL :: vtot,vtot2,vsd,vscl REAL :: wtot,wtot2,wsd,wscl REAL :: pttot,pttot2,ptsd,ptscl REAL :: qvtot,qvtot2,qvsd,qvscl REAL :: ptot,ptot2,psd,pscl REAL :: totalpoint REAL :: ampu,ampv,ampw,amppt,ampqv,ampp,ampke,enorm,escl REAL :: rateu,ratev,ratew,ratept,rateqv,ratep,rateke INTEGER :: nprocx_in, nprocy_in NAMELIST /message_passing/nproc_x, nproc_y, readsplit, nprocx_in, nprocy_in NAMELIST /prtbpara/ iensopt,inibred,iseed,iorder,isread, & soilinfl,soilfmt,iread, & uvar,vvar,wvar,ptvar,qvvar NAMELIST /input_data/ hinfmt,grdbasfn,filename NAMELIST /outpt_data/ dmp_out_joined,hdmpfmt,runnmin,grdout,basout, & varout,mstout,iceout,trbout,sfcout, & rainout,snowout,filcmprs,pertout !----------------------------------------------------------------------- ! ! Variables for mpi jobs ! !----------------------------------------------------------------------- INTEGER, PARAMETER :: fzone = 3 INTEGER :: nxlg, nylg ! global domain INTEGER :: ii,jj,ia,ja CHARACTER(LEN=256) :: tmpstr ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! !----------------------------------------------------------------------- ! ! Initializations ! !----------------------------------------------------------------------- ! CALL mpinit_proc IF(myproc == 0) WRITE(6,'(/6(/5x,a)/)') & '###############################################################', & '# #', & '# Welcome to ARPSENSIC, a program that reads in history files #', & '# generated by ARPS and produces perturbation grids variables #', & '# #', & '###############################################################' soilfmt = 1 !------------------------------------------------------------------------ ! Read namelist !------------------------------------------------------------------------ IF(myproc == 0) THEN READ(5,message_passing) WRITE(6,'(a)')'Namelist message_passing was successfully read.' write(6,message_passing) END IF CALL mpupdatei(nproc_x,1) CALL mpupdatei(nproc_y,1) CALL mpupdatei(readsplit,1) IF (mp_opt == 0) THEN nproc_x = 1 nproc_y = 1 nprocx_in = 1 nprocy_in = 1 readsplit = 0 nprocs = 1 max_fopen = 1 ELSE max_fopen = nproc_x * nproc_y END IF CALL mpinit_var IF(myproc == 0) THEN READ(5,prtbpara,END=999) write(6,prtbpara) END IF CALL mpupdatei(iensopt,1) CALL mpupdatei(inibred,1) CALL mpupdatei(iseed,1) CALL mpupdater(uvar,1) CALL mpupdater(vvar,1) CALL mpupdater(wvar,1) CALL mpupdater(ptvar,1) CALL mpupdater(qvvar,1) CALL mpupdater(iorder,1) CALL mpupdatec(soilinfl,256) CALL mpupdatei(isread,1) CALL mpupdatei(iread,1) IF(myproc == 0) THEN READ(5,input_data,END=999) write(6,input_data) END IF CALL mpupdatei(hinfmt,1) CALL mpupdatec(grdbasfn,3*256) CALL mpupdatec(filename,3*256) IF(myproc == 0) THEN READ(5,outpt_data,END=999) write(6,outpt_data) END IF CALL mpupdatei(dmp_out_joined,1) CALL mpupdatei(hdmpfmt,1) CALL mpupdater(runnmin,80) CALL mpupdatei(grdout,1) CALL mpupdatei(basout,1) CALL mpupdatei(varout,1) CALL mpupdatei(mstout,1) CALL mpupdatei(iceout,1) CALL mpupdatei(trbout,1) CALL mpupdatei(rainout,1) CALL mpupdatei(sfcout,1) CALL mpupdatei(snowout,1) CALL mpupdatei(filcmprs,1) CALL mpupdatei(pertout,1) joindmp = dmp_out_joined IF (mp_opt > 0) THEN ! should moved into mpinit_var later dumpstride = max_fopen IF (joindmp > 0) dumpstride = nprocs ! join and dump END IF GO TO 1001 999 IF(myproc == 0) PRINT *, 'ERROR in reading from input file, Program Terminated' STOP 1001 CONTINUE !------------------------------------------------------------------------ ! ! Get the dimensions of the input data files ! !------------------------------------------------------------------------ lengbf(1) = LEN_TRIM(grdbasfn(1)) IF(mp_opt > 0 .AND. readsplit <= 0) THEN tmpstr = grdbasfn(1) WRITE(grdbasfn(1),'(a,a,2i2.2)') tmpstr(1:lengbf(1)),'_',loc_x,loc_y lengbf(1) = lengbf(1) + 5 END IF IF(myproc == 0) THEN PRINT *, 'Get dimension from:',trim(grdbasfn(1)) CALL get_dims_from_data(hinfmt,grdbasfn(1)(1:lengbf(1)), & nx,ny,nz,nzsoil,nstyps, ireturn) IF (mp_opt > 0 .AND. readsplit > 0) THEN ! ! Fiddle with nx/ny, which apparently are wrong. ! nx = (nx - 3) / nproc_x + 3 ny = (ny - 3) / nproc_y + 3 END IF IF (nstyps <= 0) nstyps = 1 nstyp = nstyps ! Copy to global variabl IF( ireturn /= 0 ) THEN PRINT*,'Problem occured when trying to get dimensions from data.' PRINT*,'Program stopped.' CALL arpsstop('Problem with data.',1) END IF END IF CALL mpupdatei(nx,1) CALL mpupdatei(ny,1) CALL mpupdatei(nz,1) CALL mpupdatei(nzsoil,1) CALL mpupdatei(nstyps,1) CALL mpupdatei(nstyp,1) IF(myproc == 0) & WRITE(6,'(4(a,i5))') 'nx =',nx,', ny=',ny,', nz=',nz,' nzsoil =',nzsoil IF(myproc == 0) print*,'nstyps =', nstyps ! !----------------------------------------------------------------------- ! ! Allocate the variables and initialize the them to zero ! !----------------------------------------------------------------------- ALLOCATE(x(nx),STAT=istatus) x=0 ALLOCATE(y(ny),STAT=istatus) y=0 ALLOCATE(z(nz),STAT=istatus) z=0 ALLOCATE(zp(nx,ny,nz),STAT=istatus) zp=0 ALLOCATE(zpsoil(nx,ny,nzsoil),STAT=istatus) zpsoil=0 ALLOCATE(uprt(nx,ny,nz),STAT=istatus) uprt=0 ALLOCATE(vprt(nx,ny,nz),STAT=istatus) vprt=0 ALLOCATE(wprt(nx,ny,nz),STAT=istatus) wprt=0 ALLOCATE(ptprt(nx,ny,nz),STAT=istatus) ptprt=0 ALLOCATE(pprt(nx,ny,nz),STAT=istatus) pprt=0 ALLOCATE(qvprt(nx,ny,nz),STAT=istatus) qvprt=0 ALLOCATE(qc(nx,ny,nz),STAT=istatus) qc=0 ALLOCATE(qr(nx,ny,nz),STAT=istatus) qr=0 ALLOCATE(qi(nx,ny,nz),STAT=istatus) qi=0 ALLOCATE(qs(nx,ny,nz),STAT=istatus) qs=0 ALLOCATE(qh(nx,ny,nz),STAT=istatus) qh=0 ALLOCATE(tke(nx,ny,nz),STAT=istatus) tke=0 ALLOCATE(kmh(nx,ny,nz),STAT=istatus) kmh=0 ALLOCATE(kmv(nx,ny,nz),STAT=istatus) kmv=0 ALLOCATE(ubar(nx,ny,nz),STAT=istatus) ubar=0 ALLOCATE(vbar(nx,ny,nz),STAT=istatus) vbar=0 ALLOCATE(wbar(nx,ny,nz),STAT=istatus) wbar=0 ALLOCATE(ptbar(nx,ny,nz),STAT=istatus) ptbar=0 ALLOCATE(pbar(nx,ny,nz),STAT=istatus) pbar=0 ALLOCATE(rhobar(nx,ny,nz),STAT=istatus) rhobar=0 ALLOCATE(qvbar(nx,ny,nz),STAT=istatus) qvbar=0 ALLOCATE(soiltyp(nx,ny,nstyps),STAT=istatus) soiltyp=0 ALLOCATE(stypfrct(nx,ny,nstyps),STAT=istatus) stypfrct=0 ALLOCATE(vegtyp(nx,ny),STAT=istatus) vegtyp=0 ALLOCATE(lai(nx,ny),STAT=istatus) lai=0 ALLOCATE(roufns(nx,ny),STAT=istatus) roufns=0 ALLOCATE(veg(nx,ny),STAT=istatus) veg=0 ALLOCATE(tsoil(nx,ny,nzsoil,0:nstyps),STAT=istatus) tsoil=0 ALLOCATE(qsoil(nx,ny,nzsoil,0:nstyps),STAT=istatus) qsoil=0 ALLOCATE(wetcanp(nx,ny,0:nstyps),STAT=istatus) wetcanp=0 ALLOCATE(snowdpth(nx,ny),STAT=istatus) snowdpth=0 ALLOCATE(raing(nx,ny),STAT=istatus) raing=0 ALLOCATE(rainc(nx,ny),STAT=istatus) rainc=0 ALLOCATE(vx(nx),STAT=istatus) vx=0 ALLOCATE(vy(ny),STAT=istatus) vy=0 ALLOCATE(vz(nz),STAT=istatus) vz=0 ALLOCATE(vzp(nx,ny,nz),STAT=istatus) vzp=0 ALLOCATE(vzpsoil(nx,ny,nzsoil),STAT=istatus) vzpsoil=0 ALLOCATE(vuprt(nx,ny,nz),STAT=istatus) vuprt=0 ALLOCATE(vvprt(nx,ny,nz),STAT=istatus) vvprt=0 ALLOCATE(vwprt(nx,ny,nz),STAT=istatus) vwprt=0 ALLOCATE(vptprt(nx,ny,nz),STAT=istatus) vptprt=0 ALLOCATE(vpprt(nx,ny,nz),STAT=istatus) vpprt=0 ALLOCATE(vqvprt(nx,ny,nz),STAT=istatus) vqvprt=0 ALLOCATE(vqc(nx,ny,nz),STAT=istatus) vqc=0 ALLOCATE(vqr(nx,ny,nz),STAT=istatus) vqr=0 ALLOCATE(vqi(nx,ny,nz),STAT=istatus) vqi=0 ALLOCATE(vqs(nx,ny,nz),STAT=istatus) vqs=0 ALLOCATE(vqh(nx,ny,nz),STAT=istatus) vqh=0 ALLOCATE(vtke(nx,ny,nz),STAT=istatus) vtke=0 ALLOCATE(vkmh(nx,ny,nz),STAT=istatus) vkmh=0 ALLOCATE(vkmv(nx,ny,nz),STAT=istatus) vkmv=0 ALLOCATE(vubar(nx,ny,nz),STAT=istatus) vubar=0 ALLOCATE(vvbar(nx,ny,nz),STAT=istatus) vvbar=0 ALLOCATE(vwbar(nx,ny,nz),STAT=istatus) vwbar=0 ALLOCATE(vptbar(nx,ny,nz),STAT=istatus) vptbar=0 ALLOCATE(vpbar(nx,ny,nz),STAT=istatus) vpbar=0 ALLOCATE(vrhobar(nx,ny,nz),STAT=istatus) vrhobar=0 ALLOCATE(vqvbar(nx,ny,nz),STAT=istatus) vqvbar=0 ALLOCATE(vsoiltyp(nx,ny,nstyps),STAT=istatus) vsoiltyp=0 ALLOCATE(vstypfrct(nx,ny,nstyps),STAT=istatus) vstypfrct=0 ALLOCATE(vvegtyp(nx,ny),STAT=istatus) vvegtyp=0 ALLOCATE(vlai(nx,ny),STAT=istatus) vlai=0 ALLOCATE(vroufns(nx,ny),STAT=istatus) vroufns=0 ALLOCATE(vveg(nx,ny),STAT=istatus) vveg=0 ALLOCATE(vtsoil(nx,ny,nzsoil,0:nstyps),STAT=istatus) vtsoil=0 ALLOCATE(vqsoil(nx,ny,nzsoil,0:nstyps),STAT=istatus) vqsoil=0 ALLOCATE(vwetcanp(nx,ny,0:nstyps),STAT=istatus) vwetcanp=0 ALLOCATE(vsnowdpth(nx,ny),STAT=istatus) vsnowdpth=0 ALLOCATE(vraing(nx,ny),STAT=istatus) vraing=0 ALLOCATE(vrainc(nx,ny),STAT=istatus) vrainc=0 ALLOCATE(prcrate(nx,ny,4),STAT=istatus) prcrate=0 ALLOCATE(radfrc(nx,ny,nz),STAT=istatus) radfrc=0 ALLOCATE(radsw(nx,ny),STAT=istatus) radsw=0 ALLOCATE(rnflx(nx,ny),STAT=istatus) rnflx=0 ALLOCATE(radswnet(nx,ny),STAT=istatus) radswnet=0 ALLOCATE(radlwin(nx,ny),STAT=istatus) radlwin=0 ALLOCATE(usflx(nx,ny),STAT=istatus) usflx=0 ALLOCATE(vsflx(nx,ny),STAT=istatus) vsflx=0 ALLOCATE(ptsflx(nx,ny),STAT=istatus) ptsflx=0 ALLOCATE(qvsflx(nx,ny),STAT=istatus) qvsflx=0 ALLOCATE(tem1(nx,ny,nz),STAT=istatus) tem1=0 ALLOCATE(tem2(nx,ny,nz),STAT=istatus) tem2=0 ALLOCATE(tem3(nx,ny,nz),STAT=istatus) tem3=0 ALLOCATE(tem3dsoil(nx,ny,nzsoil),STAT=istatus) tem3dsoil=0 ALLOCATE(vtem1(nx,ny,nz),STAT=istatus) vtem1=0 ALLOCATE(vtem2(nx,ny,nz),STAT=istatus) vtem2=0 ALLOCATE(vtem3(nx,ny,nz),STAT=istatus) vtem3=0 ALLOCATE(vtem3dsoil(nx,ny,nzsoil),STAT=istatus) vtem3dsoil=0 ALLOCATE(xs(nx),STAT=istatus) xs=0 ALLOCATE(ys(ny),STAT=istatus) ys=0 ALLOCATE(zps(nx,ny,nz),STAT=istatus) zps=0 ALLOCATE(lat(nx,ny),STAT=istatus) ALLOCATE(lon(nx,ny),STAT=istatus) lat=0 lon=0 ALLOCATE(vxs(nx),STAT=istatus) vxs=0 ALLOCATE(vys(ny),STAT=istatus) vys=0 ALLOCATE(vzps(nx,ny,nz),STAT=istatus) vzps=0 ALLOCATE(dxfld(nx),STAT=istatus) dxfld=0 ALLOCATE(dyfld(ny),STAT=istatus) dyfld=0 ALLOCATE(rdxfld(nx),STAT=istatus) rdxfld=0 ALLOCATE(rdyfld(ny),STAT=istatus) rdyfld=0 ! ! mpi variables ! nxlg = (nx-fzone)*nproc_x + fzone nylg = (ny-fzone)*nproc_y + fzone totalpoint = (nxlg-3)*(nylg-3)*(nz-1) IF(myproc == 0) print *, 'totalpoint =',totalpoint 101 CONTINUE ! !----------------------------------------------------------------------- ! ! Get the name of the grid/base data set. ! !----------------------------------------------------------------------- ! IF (myproc == 0) THEN DO i=1,3 WRITE(6,'(/a,i5,a)') ' For data set ',i,':' WRITE(6,'(/a,/1x,a)')' The base set name is ',trim(grdbasfn(i)) WRITE(6,'(/a,/1x,a)')' The data set name is ',trim(filename(i)) END DO WRITE(6,'(/5x,a,a)') 'The output run name is: ', trim(runnmin) WRITE(6,'(a,i5)') ' The output data format flag is: ',hdmpfmt END IF !----------------------------------------------------------------------- ! ! Read all input data arrays (a - forecast data) ! !----------------------------------------------------------------------- ! lenfil(1) = LEN_TRIM(filename(1)) IF (mp_opt > 0 .AND. readsplit <= 0) THEN tmpstr = filename(1) WRITE(filename(1),'(2a,2i2.2)') tmpstr(1:lenfil(1)),'_',loc_x,loc_y lenfil(1) = lenfil(1) +5 END IF IF (myproc == 0) THEN WRITE(6,'(/2a/)') ' Reading file: ',filename(1)(1:lenfil(1)) WRITE(6,'(/2a/)') ' Reading file: ',grdbasfn(1)(1:lengbf(1)) END IF CALL dtaread(nx,ny,nz,nzsoil,nstyps, & hinfmt,nchin,grdbasfn(1)(1:lengbf(1)),lengbf(1), & filename(1)(1:lenfil(1)),lenfil(1),time, & x,y,z,zp,zpsoil, uprt ,vprt ,wprt ,ptprt, pprt , & qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, & ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, & soiltyp,stypfrct,vegtyp,lai,roufns,veg, & tsoil,qsoil,wetcanp,snowdpth, & raing,rainc,prcrate, & radfrc,radsw,rnflx,radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & ireturn, tem1,tem2,tem3) IF (isread == 1) THEN ! currently, no MPI function CALL readsoil(nx,ny,nzsoil,nstyps,soilinfl,dx,dy,zpsoil, & mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, & zpsoilin,tsoilin,qsoilin,wcanpin,snowdin, & tsoil,qsoil,wetcanp,snowdpth,soiltyp) END IF ! IF( ireturn /= 0 ) CALL arpsstop('dtaread errors.',1) ! IF (inibred == 1) THEN ! initial bred, preparing sd(f) utot=0.0 vtot=0.0 wtot=0.0 pttot=0.0 qvtot=0.0 ptot=0.0 utot2=0.0 vtot2=0.0 wtot2=0.0 pttot2=0.0 qvtot2=0.0 ptot2=0.0 do k=1,nz-1 do i=2,nx-2 do j=2,ny-2 utot=utot+(ubar(i,j,k)+uprt(i,j,k)) vtot=vtot+(vbar(i,j,k)+vprt(i,j,k)) wtot=wtot+(wbar(i,j,k)+wprt(i,j,k)) pttot=pttot+ptprt(i,j,k) ! pttot=pttot+ptprt(i,j,k)+ptbar(i,j,k) qvtot=qvtot+(qvbar(i,j,k)+qvprt(i,j,k)) ptot=ptot+pprt(i,j,k) ! ptot=ptot+(pprt(i,j,k)+pbar(i,j,k)) utot2=utot2+(ubar(i,j,k)+uprt(i,j,k))*(ubar(i,j,k)+uprt(i,j,k)) vtot2=vtot2+(vbar(i,j,k)+vprt(i,j,k))*(vbar(i,j,k)+vprt(i,j,k)) wtot2=wtot2+(wbar(i,j,k)+wprt(i,j,k))*(wbar(i,j,k)+wprt(i,j,k)) pttot2=pttot2+(ptprt(i,j,k)*ptprt(i,j,k)) ! pttot2=pttot2+(ptbar(i,j,k)+ptprt(i,j,k))*(ptbar(i,j,k)+ptprt(i,j,k)) qvtot2=qvtot2+(qvbar(i,j,k)+qvprt(i,j,k))*(qvbar(i,j,k)+qvprt(i,j,k)) ptot2=ptot2+(pprt(i,j,k)*pprt(i,j,k)) ! ptot2=ptot2+(pbar(i,j,k)+pprt(i,j,k))*(pbar(i,j,k)+pprt(i,j,k)) enddo enddo enddo IF (mp_opt > 0) THEN CALL mpsumr(utot, 1) CALL mpsumr(vtot, 1) CALL mpsumr(wtot, 1) CALL mpsumr(pttot, 1) CALL mpsumr(qvtot, 1) CALL mpsumr(ptot, 1) CALL mpsumr(utot2, 1) CALL mpsumr(vtot2, 1) CALL mpsumr(wtot2, 1) CALL mpsumr(pttot2, 1) CALL mpsumr(qvtot2, 1) CALL mpsumr(ptot2, 1) END IF utot=utot/totalpoint vtot=utot/totalpoint wtot=wtot/totalpoint pttot=pttot/totalpoint qvtot=qvtot/totalpoint ptot=ptot/totalpoint usd=sqrt(utot2/totalpoint-utot*utot) vsd=sqrt(vtot2/totalpoint-vtot*vtot) wsd=sqrt(wtot2/totalpoint-wtot*wtot) ptsd=sqrt(pttot2/totalpoint-pttot*pttot) qvsd=sqrt(qvtot2/totalpoint-qvtot*qvtot) psd=sqrt(ptot2/totalpoint-ptot*ptot) IF (myproc == 0) THEN print *,'totalpoint=',totalpoint print *,'variance info:' print *,'usd=',usd,' utot=',utot,' utot2=',utot2/totalpoint print *,'vsd=',vsd,' vtot=',vtot,' vtot2=',vtot2/totalpoint print *,'wsd=',wsd,' wtot=',wtot,' wtot2=',wtot2/totalpoint print *,'ptsd=',ptsd,' pttot=',pttot,' pttot2=',pttot2/totalpoint print *,'qvsd=',qvsd,' qvtot=',qvtot,' qvtot2=',qvtot2/totalpoint print *,'psd=',psd,' ptot=',ptot,' ptot2=',ptot2/totalpoint write(12,*) 'utot2,vtot2,wtot2,pttot2,qvtot2,ptot2,enorm - total mean' write(12,'(7e11.4)') sqrt(utot2/totalpoint),sqrt(vtot2/totalpoint), & sqrt(wtot2/totalpoint),sqrt(pttot2/totalpoint), & sqrt(qvtot2/totalpoint),sqrt(ptot2/totalpoint), & sqrt(0.5*(utot2+vtot2+wtot2)/totalpoint) endif ! set initial (random) amplitude for each variable ampu = uvar ! provided from input file ampv = vvar ampw = wvar amppt= ptvar ampqv= qvvar ampp = 0.0 IF (myproc == 0) THEN print *,'Initial amplitudes (for random perturbation):' print *,'ampu =',ampu print *,'ampv =',ampv print *,'ampw =',ampw print *,'amppt=',amppt print *,'ampqv=',ampqv print *,'ampp=',ampp endif ELSE ampu = uvar ! provided from input file ampv = vvar ampw = wvar amppt= ptvar ampqv= qvvar ampp = 0.0 IF (myproc == 0) THEN print *,'Perturbation switches ( >0 on; =0 off):' print *,'pert_u =',ampu print *,'pert_v =',ampv print *,'pert_w =',ampw print *,'pert_pt=',amppt print *,'pert_qv=',ampqv print *,'pert_p =',ampp endif END IF curtim=time fcrnam=runname ifproj=mapproj fscale=sclfct flatnot(1)=trulat1 flatnot(2)=trulat2 ftrulon=trulon fdx=x(3)-x(2) fdy=y(3)-y(2) fctrlat=ctrlat fctrlon=ctrlon CALL setmapr(ifproj,fscale,flatnot,ftrulon) CALL lltoxy(1,1,fctrlat,fctrlon,xctr,yctr) fx0=xctr-fdx*((nxlg-3)/2) fy0=yctr-fdy*((nylg-3)/2) CALL setorig(1,fx0,fy0) ! !----------------------------------------------------------------------- ! ! Get verification data (b - analysis (or gridded obs.) data) ! !----------------------------------------------------------------------- ! ! Set the gridread parameter to 0 so that the verification ! grid/base file will be read. ! !----------------------------------------------------------------------- ! CALL setgbrd (0) ! !----------------------------------------------------------------------- ! ! Read in the verification data. ! !----------------------------------------------------------------------- ! lenfil(2) = len_trim(filename(2)) lengbf(2) = len_trim(grdbasfn(2)) IF (mp_opt > 0 .AND. readsplit <= 0) THEN tmpstr = filename(2) WRITE(filename(2),'(2a,2i2.2)') tmpstr(1:lenfil(2)),'_',loc_x,loc_y lenfil(2) = lenfil(2) + 5 tmpstr = grdbasfn(2) WRITE(grdbasfn(2),'(2a,2i2.2)') tmpstr(1:lengbf(2)),'_',loc_x,loc_y lengbf(2) = lengbf(2) + 5 END IF IF (myproc == 0) THEN WRITE(6,'(/2a/)') ' Reading file: ',filename(2)(1:lenfil(2)) WRITE(6,'(/2a/)') ' Reading file: ',grdbasfn(2)(1:lengbf(2)) END IF CALL dtaread(nx,ny,nz,nzsoil,nstyps, & hinfmt,nchin,grdbasfn(2)(1:lengbf(2)),lengbf(2), & filename(2)(1:lenfil(2)),lenfil(2),time, & vx,vy,vz,vzp,vzpsoil,vuprt,vvprt,vwprt,vptprt,vpprt, & vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv, & vubar, vvbar, vwbar, vptbar, vpbar, vrhobar, vqvbar, & vsoiltyp,vstypfrct,vvegtyp,vlai,vroufns,vveg, & vtsoil,vqsoil,vwetcanp,vsnowdpth, & raing,rainc,prcrate, & radfrc,radsw,rnflx,radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & ireturn, vtem1,vtem2,vtem3) IF (isread == 1) THEN CALL readsoil(nx,ny,nzsoil,nstyps,soilinfl,dx,dy,vzpsoil, & mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, & zpsoilin,tsoilin,qsoilin,wcanpin,snowdin, & vtsoil,vqsoil,vwetcanp,vsnowdpth,vsoiltyp) END IF IF( ireturn /= 0 ) CALL arpsstop('dtaread errors.',1) ! IF (vnx /= nx.OR.vny /= ny.OR.vnz /= nz) THEN ! PRINT *,'nx,ny,nz,','=/=','vnx,vny,vnz' ! PRINT *,nx,ny,nz,'=/=',vnx,vny,vnz ! PRINT *, 'forced to stop' ! STOP ! END IF ! don't need ! ivproj=mapproj ! vscale=sclfct ! vlatnot(1)=trulat1 ! vlatnot(2)=trulat2 ! vtrulon=trulon ! vdx=vx(3)-vx(2) ! vdy=vy(3)-vy(2) ! vctrlat=ctrlat ! vctrlon=ctrlon ! CALL setmapr(ivproj,vscale,vlatnot,vtrulon) ! CALL lltoxy(1,1,vctrlat,vctrlon,xctr,yctr) ! vx0=xctr-vdx*((vnx-3)/2) ! vy0=yctr-vdy*((vny-3)/2) ! CALL setorig(1,vx0,vy0) ! ! IF (fx0 == vx0.AND.fy0 == vy0.AND. & ! flatnot(1) == vlatnot(1).AND.flatnot(2) == vlatnot(2).AND. & ! ftrulon == vtrulon.AND.ifproj == ivproj .AND. & ! fscale == vscale ) THEN ! IF(myproc == 0) PRINT *, 'Grids 1 and 2 shares a common coordinate system' ! ELSE ! IF(myproc == 0) THEN ! PRINT *, 'Grids 1/2 are different, CHECK the PROGRAM or data' ! PRINT *, 'Forced to STOP' ! ENDIF ! STOP ! END IF ! IF (inibred == 0) THEN ! !----------------------------------------------------------------------- ! ! Find difference = forecast - verification ! a=a-b (note the forth parameter is 0) ! To reduce memory requirements, the difference fields are ! written to the same arrays as the interpolated fields. ! !----------------------------------------------------------------------- ! CALL prtfield(nx,ny,nz,nzsoil,0, & uprt, vprt, wprt, ptprt, pprt, & qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, & ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, & tsoil,qsoil,wetcanp, & raing,rainc, & vuprt, vvprt, vwprt, vptprt, vpprt, & vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv, & vubar, vvbar, vwbar, vptbar, vpbar, vrhobar, vqvbar, & vtsoil,vqsoil,vwetcanp, & vraing,vrainc, & uprt, vprt, wprt, ptprt, pprt, & qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, & tsoil,qsoil,wetcanp, & raing,rainc, & tem1,tem3dsoil,ireturn ) ! write out 3D perturbation fields for plotting purpose IF(pertout > 0) THEN DO i=0,nprocs-1,max_fopen IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN filnam = filename(1)(1:lenfil(1)-9)//'ptpert'// & filename(1)(lenfil(1)-5:lenfil(1)) ! Here, further action is needed for split file name (KONG) q3dname='pt_pert' q3dunit='K' CALL bindump3d(nx,ny,nz,trim(filnam),q3dname,q3dunit,ptprt, & 1,1,1) filnam = filename(1)(1:lenfil(1)-9)//'qvpert'// & filename(1)(lenfil(1)-5:lenfil(1)) q3dname='qv_pert' q3dunit='g/kg' CALL bindump3d(nx,ny,nz,trim(filnam),q3dname,q3dunit,1e3*qvprt, & 1,1,1) filnam = filename(1)(1:lenfil(1)-9)//'u_pert'// & filename(1)(lenfil(1)-5:lenfil(1)) q3dname='u_pert' q3dunit='m/s' CALL bindump3d(nx,ny,nz,trim(filnam),q3dname,q3dunit,uprt, & 0,1,1) filnam = filename(1)(1:lenfil(1)-9)//'v_pert'// & filename(1)(lenfil(1)-5:lenfil(1)) q3dname='v_pert' q3dunit='m/s' CALL bindump3d(nx,ny,nz,trim(filnam),q3dname,q3dunit,vprt, & 1,0,1) filnam = filename(1)(1:lenfil(1)-9)//'w_pert'// & filename(1)(lenfil(1)-5:lenfil(1)) q3dname='w_pe' CALL bindump3d(nx,ny,nz,trim(filnam),q3dname,q3dunit,wprt, & 1,1,0) filnam = filename(1)(1:lenfil(1)-9)//'p_pert'// & filename(1)(lenfil(1)-5:lenfil(1)) q3dname='p_pert' q3dunit='Pa' CALL bindump3d(nx,ny,nz,trim(filnam),q3dname,q3dunit,pprt, & 1,1,1) END IF IF (mp_opt > 0) CALL mpbarrier END DO END IF ! End writing 3D perturbation fields ! ! Select perturbation fields based on switches IF(ampu < 1e-5) uprt = 0.0 IF(ampv < 1e-5) vprt = 0.0 IF(ampw < 1e-5) wprt = 0.0 IF(amppt< 1e-5) ptprt= 0.0 IF(ampqv< 1e-5) qvprt= 0.0 IF(ampp < 1e-5) pprt = 0.0 utot2=0.0 vtot2=0.0 wtot2=0.0 pttot2=0.0 qvtot2=0.0 ptot2=0.0 do k=1,nz-1 do i=2,nx-2 do j=2,ny-2 utot2=utot2+uprt(i,j,k)*uprt(i,j,k) vtot2=vtot2+vprt(i,j,k)*vprt(i,j,k) wtot2=wtot2+wprt(i,j,k)*wprt(i,j,k) pttot2=pttot2+ptprt(i,j,k)*ptprt(i,j,k) qvtot2=qvtot2+qvprt(i,j,k)*qvprt(i,j,k) ptot2=ptot2+pprt(i,j,k)*pprt(i,j,k) enddo enddo enddo IF (mp_opt > 0) THEN CALL mpsumr(utot2, 1) CALL mpsumr(vtot2, 1) CALL mpsumr(wtot2, 1) CALL mpsumr(pttot2, 1) CALL mpsumr(qvtot2, 1) CALL mpsumr(ptot2, 1) END IF usd=sqrt(utot2/totalpoint) vsd=sqrt(vtot2/totalpoint) wsd=sqrt(wtot2/totalpoint) ptsd=sqrt(pttot2/totalpoint) qvsd=sqrt(qvtot2/totalpoint) psd=sqrt(ptot2/totalpoint) enorm=sqrt(0.5*(utot2+vtot2+wtot2)/totalpoint) IF (myproc == 0) THEN print *,'Perturbation norms before scaling:' print *,'unorm =',usd print *,'vnorm =',vsd print *,'wnorm =',wsd print *,'ptnorm=',ptsd print *,'qvnorm=',qvsd print *,'Pnorm =',psd print *,'KEnorm =',enorm write(13,*) 'usd,vsd,wsd,ptsd,qvsd,psd,enorm' write(13,'(7e11.4)') usd,vsd,wsd,ptsd,qvsd,psd,enorm endif IF (iensopt == 1) THEN ! scaling the perturbation in breeding IC IF (myproc == 0) THEN read(10,*) ampu,ampv,ampw,amppt,ampqv,ampp,ampke ! now hold initial norms endif CALL mpupdater(ampu,1) CALL mpupdater(ampv,1) CALL mpupdater(ampw,1) CALL mpupdater(amppt,1) CALL mpupdater(ampqv,1) CALL mpupdater(ampp,1) CALL mpupdater(ampke,1) ! for rescaling ! or assign valuess here manually for rescaling IF (myproc == 0) THEN print *,'Initial norms read back:' print *,ampu,ampv,ampw,amppt,ampqv,ampp,ampke endif if(usd > 0.0) uscl = ampu/usd if(vsd > 0.0) vscl = ampv/vsd if(wsd > 0.0) wscl = ampw/wsd if(ptsd > 0.0) ptscl = amppt/ptsd if(qvsd > 0.0) qvscl = ampqv/qvsd if(psd > 0.0) pscl = ampp/psd if(enorm > 0.0) escl=ampke/enorm rateu=-1.0 ! special value (no meaning) ratev=-1.0 ratew=-1.0 ratept=-1.0 rateqv=-1.0 ratep=-1.0 rateke=-1.0 if(uscl > 0.0) rateu=1.0/uscl if(vscl > 0.0) ratev=1.0/vscl if(wscl > 0.0) ratew=1.0/wscl if(ptscl > 0.0) ratept=1.0/ptscl if(qvscl > 0.0) rateqv=1.0/qvscl if(pscl > 0.0) ratep=1.0/pscl if(escl > 0.0) rateke=1.0/escl IF(myproc == 0) THEN write(11,'(7f11.4)') usd,vsd,wsd,ptsd,qvsd,psd,enorm write(11,'(7f11.4)') uscl,vscl,wscl,ptscl,qvscl,pscl,escl write(11,'(7f11.4)') rateu,ratev,ratew,ratept,rateqv,ratep,rateke print *,'totalpoint=',totalpoint print *,'usd=',usd,' uscl=',uscl print *,'vsd=',vsd,' vscl=',vscl print *,'wsd=',wsd,' wscl=',wscl print *,'ptsd=',ptsd,' ptscl=',ptscl print *,'qvsd=',qvsd,' qvscl=',qvscl print *,'psd=',psd,' pscl=',pscl endif uprt = uscl*uprt vprt = vscl*vprt wprt = wscl*wprt ptprt= ptscl*ptprt qvprt= qvscl*qvprt pprt = pscl*pprt ELSE IF(myproc == 0) THEN print *,'Perturbation norms after scaling:' print *,'unorm =',abs(iorder)*usd print *,'vnorm =',abs(iorder)*vsd print *,'wnorm =',abs(iorder)*wsd print *,'ptnorm=',abs(iorder)*ptsd print *,'qvnorm=',abs(iorder)*qvsd print *,'Pnorm =',abs(iorder)*psd print *,'KEnorm =',abs(iorder)*enorm write(14,*) 'usd,vsd,wsd,ptsd,qvsd,psd,enorm - after scaling' write(14,'(7e11.4)') abs(iorder)*usd,abs(iorder)*vsd,abs(iorder)*wsd, & abs(iorder)*ptsd,abs(iorder)*qvsd,abs(iorder)*psd,abs(iorder)*enorm endif END IF ELSE IF (inibred ==1) THEN ! generate random perturbations uprt=0.0 vprt=0.0 wprt=0.0 ptprt=0.0 qvprt=0.0 pprt=0.0 if(iseed < 0) then ! computer select iseed CALL DATE_AND_TIME(VALUES=idate) iseed=idate(8) IF(myproc == 0) print *,'Machine decided ISEED=',iseed IF(myproc == 0) write(98,*) iseed else if (iseed == 9999) then IF(myproc == 0) read(98,*) iseed IF(myproc == 0) print *,'Read in ISEED: ',iseed CALL mpupdatei(iseed,1) endif IF(ampu > 1e-10) THEN IF(myproc == 0) print *,'U-PERT initial iseed =',iseed call normalrand(iseed,nx*ny*nz,tem1) do k=1,nz-1 do i=1,nx do j=1,ny-1 uprt(i,j,k) = ampu*tem1(i,j,k) enddo enddo enddo END IF IF(ampv > 1e-10) THEN IF(myproc == 0) print *,'V-PERT initial iseed =',iseed call normalrand(iseed,nx*ny*nz,tem1) do k=1,nz-1 do i=1,nx-1 do j=1,ny vprt(i,j,k) = ampv*tem1(i,j,k) enddo enddo enddo END IF IF(ampw > 1e-10) THEN IF(myproc == 0) print *,'W-PERT initial iseed =',iseed call normalrand(iseed,nx*ny*nz,tem1) do k=1,nz-1 do i=1,nx-1 do j=1,ny-1 wprt(i,j,k) = ampw*tem1(i,j,k) enddo enddo enddo END IF IF(amppt > 1e-10) THEN IF(myproc == 0) print *,'PT-PERT initial iseed =',iseed call normalrand(iseed,nx*ny*nz,tem1) do k=1,nz-1 do i=1,nx-1 do j=1,ny-1 ptprt(i,j,k) = amppt*tem1(i,j,k) enddo enddo enddo END IF IF(ampqv > 1e-10) THEN IF(myproc == 0) print *,'QV-PERT initial iseed =',iseed call normalrand(iseed,nx*ny*nz,tem1) do k=1,nz-1 do i=1,nx-1 do j=1,ny-1 qvprt(i,j,k) = ampqv*tem1(i,j,k) enddo enddo enddo END IF IF(ampp > 1e-10) THEN IF(myproc == 0) print *,'P-PERT initial iseed =',iseed call normalrand(iseed,nx*ny*nz,tem1) do k=1,nz-1 do i=1,nx-1 do j=1,ny-1 pprt(i,j,k) = ampp*tem1(i,j,k) enddo enddo enddo END IF tem1=0.0 qc=0.0 qr=0.0 qi=0.0 qs=0.0 qh=0.0 tke=0.0 kmh=0.0 kmv=0.0 tsoil=0.0 qsoil=0.0 wetcanp=0.0 raing=0.0 rainc=0.0 ! Calculate and save initial perturbations (norm) utot2=0.0 vtot2=0.0 wtot2=0.0 pttot2=0.0 qvtot2=0.0 ptot2=0.0 do k=1,nz-1 do i=2,nx-2 do j=2,ny-2 utot2=utot2+uprt(i,j,k)*uprt(i,j,k) vtot2=vtot2+vprt(i,j,k)*vprt(i,j,k) wtot2=wtot2+wprt(i,j,k)*wprt(i,j,k) pttot2=pttot2+ptprt(i,j,k)*ptprt(i,j,k) qvtot2=qvtot2+qvprt(i,j,k)*qvprt(i,j,k) ptot2=ptot2+pprt(i,j,k)*pprt(i,j,k) enddo enddo enddo IF (mp_opt > 0) THEN CALL mpsumr(utot2, 1) CALL mpsumr(vtot2, 1) CALL mpsumr(wtot2, 1) CALL mpsumr(pttot2, 1) CALL mpsumr(qvtot2, 1) CALL mpsumr(ptot2, 1) END IF usd=sqrt(utot2/totalpoint) vsd=sqrt(vtot2/totalpoint) wsd=sqrt(wtot2/totalpoint) ptsd=sqrt(pttot2/totalpoint) qvsd=sqrt(qvtot2/totalpoint) psd=sqrt(ptot2/totalpoint) enorm=sqrt(0.5*(utot2+vtot2+wtot2)/totalpoint) IF(myproc == 0) THEN print *,'Initial (random) perturbation norms:' print *,'unorm =',usd print *,'vnorm =',vsd print *,'wnorm =',wsd print *,'ptnorm=',ptsd print *,'qvnorm=',qvsd print *,'Pnorm =',psd print *,'KEnorm =',enorm write(10,'(7f11.4)') usd,vsd,wsd,ptsd,qvsd,psd,enorm endif END IF ! !----------------------------------------------------------------------- ! ! Get the name of the CONTROL data set, field c or c=b (if iread=0) ! (note: It is needed only if control .ne. verification) ! !----------------------------------------------------------------------- ! IF (iread /= 0) THEN ! !----------------------------------------------------------------------- ! ! Read all input data arrays ! !----------------------------------------------------------------------- ! CALL setgbrd (0) lenfil(3) = len_trim(filename(3)) lengbf(3) = len_trim(grdbasfn(3)) IF (mp_opt > 0 .AND. readsplit <= 0) THEN tmpstr = filename(3) WRITE(filename(3),'(2a,2i2.2)') tmpstr(1:lenfil(3)),'_',loc_x,loc_y lenfil(3) = lenfil(3) + 5 tmpstr = grdbasfn(3) WRITE(grdbasfn(3),'(2a,2i2.2)') tmpstr(1:lengbf(3)),'_',loc_x,loc_y lengbf(3) = lengbf(3) + 5 END IF IF (myproc == 0) THEN WRITE(6,'(/2a/)') ' Reading file: ',filename(3)(1:lenfil(3)) WRITE(6,'(/2a/)') ' Reading file: ',grdbasfn(3)(1:lengbf(3)) END IF CALL dtaread(nx,ny,nz,nzsoil,nstyps, & hinfmt,nchin,grdbasfn(3)(1:lengbf(3)),lengbf(3), & filename(3)(1:lenfil(3)),lenfil(3),time, & vx,vy,vz,vzp,vzpsoil,vuprt,vvprt,vwprt,vptprt,vpprt, & vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv, & vubar, vvbar, vwbar, vptbar, vpbar, vrhobar, vqvbar, & vsoiltyp,vstypfrct,vvegtyp,vlai,vroufns,vveg, & vtsoil,vqsoil,vwetcanp, vsnowdpth, & raing,rainc,prcrate, & radfrc,radsw,rnflx,radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & ireturn, vtem1,vtem2,vtem3) IF (isread == 1) THEN CALL readsoil(nx,ny,nzsoil,nstyps,soilinfl,dx,dy,vzpsoil, & mapproj,trulat1,trulat2,trulon,sclfct,ctrlat,ctrlon, & zpsoilin,tsoilin,qsoilin,wcanpin,snowdin, & vtsoil, vqsoil, vwetcanp,vsnowdpth,vsoiltyp) END IF IF( ireturn /= 0 ) CALL arpsstop('dtaread errors.',1) ! IF (vnx /= nx.OR.vny /= ny.OR.vnz /= nz) THEN ! PRINT *,'nx,ny,nz','=/=','vnx,vny,vnz' ! PRINT *,nx,ny,nz,'=/=',vnx,vny,vnz ! PRINT *, 'forced to stop' ! STOP ! END IF ! ivproj=mapproj ! vscale=sclfct ! vlatnot(1)=trulat1 ! vlatnot(2)=trulat2 ! vtrulon=trulon ! vdx=vx(3)-vx(2) ! vdy=vy(3)-vy(2) ! vctrlat=ctrlat ! vctrlon=ctrlon ! CALL setmapr(ivproj,vscale,vlatnot,vtrulon) ! CALL lltoxy(1,1,vctrlat,vctrlon,xctr,yctr) ! vx0=xctr-vdx*((vnx-3)/2) ! vy0=yctr-vdy*((vny-3)/2) ! CALL setorig(1,vx0,vy0) !! !! IF (abs(fx0-vx0) <=1.0 .AND.abs(fy0-vy0) <=1.0 .AND. & ! IF (fx0 == vx0.AND.fy0 == vy0.AND. & ! flatnot(1) == vlatnot(1).AND.flatnot(2) == vlatnot(2).AND. & ! ftrulon == vtrulon.AND.ifproj == ivproj .AND. & ! fscale == vscale ) THEN ! PRINT *, 'Grids 2 and 3 shares a common coordinate system' ! ELSE ! PRINT *, 'Grids 2/3 are different, CHECK the PROFRAM or data' ! PRINT *, 'Forced to STOP' ! STOP ! END IF END IF !----------------------------------------------------------------------- ! ! Set output variables to forecast coordinates ! !----------------------------------------------------------------------- ! curtim=time ! mapproj=ivproj ! sclfct=vscale ! trulat1=vlatnot(1) ! trulat2=vlatnot(2) ! trulon=vtrulon ! ctrlat=vctrlat ! ctrlon=vctrlon ! !----------------------------------------------------------------------- ! ! Find ptb field = filed c + a/n (note iorder=/=0) ! (b=b+a/n in the code) ! To reduce memory requirements, the resulted fields are ! written to the same arrays as the control fields. ! !----------------------------------------------------------------------- ! IF (abs(iorder) > 1e-5) THEN CALL prtfield(nx,ny,nz,nzsoil,iorder, & vuprt, vvprt, vwprt, vptprt, vpprt, & vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv, & vubar, vvbar, vwbar, vptbar, vpbar, vrhobar, vqvbar, & vtsoil,vqsoil,vwetcanp, & vraing,vrainc, & uprt, vprt, wprt, ptprt, pprt, & qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, & ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, & tsoil,qsoil,wetcanp, & raing,rainc, & vuprt, vvprt, vwprt, vptprt, vpprt, & vqvprt, vqc, vqr, vqi, vqs, vqh, vtke,vkmh,vkmv, & vtsoil,vqsoil,vwetcanp, & vraing,vrainc, & vtem1,vtem3dsoil,ireturn ) END IF !----------------------------------------------------------------------- ! Assign the average of soil variables to every soil type !----------------------------------------------------------------------- CALL dhslini(nx,ny,nzsoil,nstyps, & tsoil,qsoil,wetcanp) ! !----------------------------------------------------------------------- ! ! Get output info ! !----------------------------------------------------------------------- ! runname=runnmin ! !----------------------------------------------------------------------- ! ! Find out the number of characters to be used to construct file ! names. ! !----------------------------------------------------------------------- ! CALL gtlfnkey( runname, lfnkey ) ! !----------------------------------------------------------------------- ! ! Find out the number of characters to be used to construct file ! names. ! !----------------------------------------------------------------------- ! CALL gtlfnkey( runname, lfnkey ) ! IF (hdmpfmt == 10.AND.grbpkbit == 0) THEN grbpkbit=16 END IF ! !----------------------------------------------------------------------- ! ! Set control parameters for ! grid, base state, moisture, and ice variable dumping. ! !----------------------------------------------------------------------- ! varout=1 totout = totin grdout = grdin basout = basin varout = varin mstout = mstin rainout = rainin prcout = prcin iceout = icein tkeout = tkein trbout = trbin sfcout = sfcin snowout = snowin landout = landin radout = radin flxout = flxin CALL gtbasfn(runname(1:lfnkey),'./',2,hdmpfmt,mgrid,nestgrd, & grdbasfn(1), lengbf(1)) IF(myproc == 0) & WRITE(6,'(/1x,a,a)') & 'Output grid/base state file is ', grdbasfn(1)(1:lengbf(1)) nchdmp = 80 grdbas = 1 ! Dump out grd and base state arrays only DO k=1,nz DO j=1,ny DO i=1,nx vuprt(i,j,k)=vubar(i,j,k)+vuprt(i,j,k) vvprt(i,j,k)=vvbar(i,j,k)+vvprt(i,j,k) vwprt(i,j,k)=vwbar(i,j,k)+vwprt(i,j,k) vqvprt(i,j,k)=vqvbar(i,j,k)+vqvprt(i,j,k) END DO END DO END DO IF (abs(iorder) > 1e-5) THEN DO k=1,nz DO j=1,ny DO i=1,nx IF (k /= 1) THEN IF ((vpprt(i,j,k)+vpbar(i,j,k)) >= & (vpprt(i,j,k-1)+vpbar(i,j,k-1)) ) THEN vpprt(i,j,k)=vpprt(i,j,k-1)+vpbar(i,j,k-1)-vpbar(i,j,k) & -(vpbar(i,j,k-1)-vpbar(i,j,k))*0.001 END IF END IF IF (vqvprt(i,j,k) <= 1.0E-8) THEN vqvprt(i,j,k)=1.0E-8 END IF END DO END DO END DO END IF IF (myproc == 0) & PRINT*,'Writing base state history dump ', grdbasfn(1)(1:lengbf(1)) ! blocking inserted for ordering i/o for message passing DO i=0,nprocs-1,max_fopen IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN CALL dtadump(nx,ny,nz,nzsoil,nstyps, & hdmpfmt,nchdmp,grdbasfn(1)(1:lengbf(1)),grdbas,filcmprs, & vuprt,vvprt,vwprt,vptprt,vpprt, & vqvprt,vqc,vqr,vqi,vqs,vqh,vtke,vkmh,vkmv, & vubar,vvbar,vwbar,vptbar,vpbar,vrhobar,vqvbar, & vx,vy,vz,vzp,vzpsoil, & vsoiltyp,vstypfrct,vvegtyp,vlai,vroufns,vveg, & vtsoil,vqsoil,vwetcanp,vsnowdpth, & vraing,vrainc,prcrate, & radfrc,radsw,rnflx,radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & vtem1,vtem2,vtem3) END IF IF (mp_opt > 0) CALL mpbarrier END DO ! !----------------------------------------------------------------------- ! ! Find a unique name hdmpfn(1:ldmpf) for history dump data set ! at time 'curtim'. ! !----------------------------------------------------------------------- ! CALL gtdmpfn(runname(1:lfnkey),'./',2, & curtim,hdmpfmt, & mgrid,nestgrd, hdmpfn, ldmpf) IF(myproc == 0) & WRITE(6,'(/1x,a,f10.0,a,a)') & 'Output file at time ',curtim,' (s) is ', hdmpfn(1:ldmpf) grdbas = 0 ! Not just dump out grd and base state arrays only ! blocking inserted for ordering i/o for message passing DO i=0,nprocs-1,max_fopen IF(myproc >= i.AND.myproc <= i+max_fopen-1)THEN CALL dtadump(nx,ny,nz,nzsoil,nstyps, & hdmpfmt,nchdmp,hdmpfn(1:ldmpf),grdbas,filcmprs, & vuprt,vvprt,vwprt,vptprt,vpprt, & vqvprt,vqc,vqr,vqi,vqs,vqh,vtke,vkmh,vkmv, & vubar,vvbar,vwbar,vptbar,vpbar,vrhobar,vqvbar, & vx,vy,vz,vzp,vzpsoil, & vsoiltyp,vstypfrct,vvegtyp,vlai,vroufns,vveg, & vtsoil,qsoil,vwetcanp,vsnowdpth, & vraing,vrainc,prcrate, & radfrc,radsw,rnflx,radswnet,radlwin, & usflx,vsflx,ptsflx,qvsflx, & vtem1,vtem2,vtem3) END IF IF (mp_opt > 0) CALL mpbarrier END DO STOP END PROGRAM arpsensic ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE PRTFIELD ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE prtfield(nx,ny,nz,nzsoil,nscale, & 2,6 uprt, vprt, wprt, ptprt, pprt, & qvprt, qc, qr, qi, qs, qh, tke,kmh,kmv, & ubar, vbar, wbar, ptbar, pbar, rhobar, qvbar, & tsoil,qsoil,wetcanp, & raing,rainc, & auprt, avprt, awprt, aptprt, apprt, & aqvprt, aqc, aqr, aqi, aqs, aqh, atke,akmh,akmv, & aubar, avbar, awbar, aptbar, apbar, arhobar, aqvbar, & atsoil,aqsoil,awetcanp, & araing,arainc, & duprt, dvprt, dwprt, dptprt, dpprt, & dqvprt, dqc, dqr, dqi, dqs, dqh, dtke,dkmh,dkmv, & dtsoil,dqsoil,dwetcanp, & draing,drainc, & tem1,tem3dsoil,ireturn) ! ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! (1) nscale=0, Subtract the forecast fields from the verification ! fields (names beginning with "a") and output to the difference ! fields (names beginning with "d"). d=( )bar+( )-abar-a ! (2) nscale=/=0. Generate perturbation fields d=( )+a*nscale with ! the bar arraies being neglected. ! The input ( ( ) ) and output ( d ) ! fields may share the same storage location. For this subroutine ! it is assumed the forecast and corresponding verification ! data are at the same physical location, however, the physical ! location may differ between variables. That is uprt and auprt ! are at the same location, but that may differ from pprt and apprt. ! !----------------------------------------------------------------------- ! ! AUTHOR: Keith Brewster Ou School of Meteorology. April 1992 ! ! MODIFICATION HISTORY: ! 14 May 1992 (KB) changed from arps2.5 to arps3.0 ! 03 Aug 1992 (KB) updated to account for changes in arps3.0 ! ! 09/07/1995 (KB) ! Added differencing of surface (soil) fields. ! ! 15/05/1998 (DH) ! converted from the difference scheme to the current multi-purpose ! version used in ARPSENS group. ! ! 15/12/1998 (DH) ! Added the 3-Dimensionity of surface (soil) fields. ! ! 05/28/2002 (J. Brotzge) ! Added tsoil/qsoil to accomodate new soil schemes. ! !----------------------------------------------------------------------- ! ! INPUT : ! nx,ny,nz,nzsoil Array dimensions for forecast field. ! ! FORECAST FIELDS: ! ! uprt perturbation x component of velocity (m/s) ! vprt perturbation y component of velocity (m/s) ! wprt perturbation vertical component of velocity in Cartesian ! coordinates (m/s). ! ! ptprt perturbation potential temperature (K) ! pprt perturbation pressure (Pascal) ! ! qvprt perturbation 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) ! tke Turbulent Kinetic Energy ((m/s)**2) ! kmh Horizontal turbulent mixing coefficient (m**2/s) ! kmv Vertical turbulent mixing coefficient (m**2/s) ! ! tsoil Soil temperature (K) ! qsoil Soil moisture (m**3/m**3) ! wetcanp Canopy water amount ! ! raing Grid supersaturation rain ! rainc Cumulus convective rain ! ! INTERPOLATED VERIFICATION FIELDS: ! ! auprt perturbation x component of velocity (m/s) ! avprt perturbation y component of velocity (m/s) ! awprt perturbation vertical component of velocity in Cartesian ! coordinates (m/s). ! ! aptprt perturbation potential temperature (K) ! apprt perturbation pressure (Pascal) ! ! aqvprt perturbation water vapor mixing ratio (kg/kg) ! aqc Cloud water mixing ratio (kg/kg) ! aqr Rainwater mixing ratio (kg/kg) ! aqi Cloud ice mixing ratio (kg/kg) ! aqs Snow mixing ratio (kg/kg) ! aqh Hail mixing ratio (kg/kg) ! ! aubar Base state x velocity component (m/s) ! avbar Base state y velocity component (m/s) ! awbar Base state z velocity component (m/s) ! aptbar Base state potential temperature (K) ! apbar Base state pressure (Pascal) ! arhobar Base state density (kg/m**3) ! aqvbar Base state water vapor mixing ratio (kg/kg) ! ! atsoil Soil temperature (K) ! aqsoil Soil moisture (m**3/m**3) ! awetcanp Canopy water amount ! ! araing Grid supersaturation rain ! arainc Cumulus convective rain ! ! OUTPUT : ! ! DIFFERENCE FIELDS (may share storage with forecast fields ! or interpolated fields in calling program): ! ! duprt perturbation x component of velocity (m/s) ! dvprt perturbation y component of velocity (m/s) ! dwprt perturbation vertical component of velocity in Cartesian ! coordinates (m/s). ! ! dptprt perturbation potential temperature (K) ! dpprt perturbation pressure (Pascal) ! ! dqvprt perturbation water vapor mixing ratio (kg/kg) ! dqc Cloud water mixing ratio (kg/kg) ! dqr Rainwater mixing ratio (kg/kg) ! dqi Cloud ice mixing ratio (kg/kg) ! dqs Snow mixing ratio (kg/kg) ! dqh Hail mixing ratio (kg/kg) ! ! dtsoil Soil temperature (K) ! dqsoil Soil moisture (m**3/m**3) ! dwetcanp Canopy water amount ! ! draing Grid supersaturation rain ! drainc Cumulus convective rain ! ! tem1 Work array ! tem3dsoil Work array ! !----------------------------------------------------------------------- ! ! Variable Declarations: ! !----------------------------------------------------------------------- ! IMPLICIT NONE INCLUDE 'mp.inc' INTEGER :: nx,ny,nz,nzsoil ! 4 dimensions of array REAL :: nscale ! !----------------------------------------------------------------------- ! ! Model Arrays ! !----------------------------------------------------------------------- ! REAL :: uprt (nx,ny,nz) ! Perturbation u-velocity (m/s) REAL :: vprt (nx,ny,nz) ! Perturbation v-velocity (m/s) REAL :: wprt (nx,ny,nz) ! Perturbation w-velocity (m/s) REAL :: ptprt (nx,ny,nz) ! Perturbation potential temperature (K) REAL :: pprt (nx,ny,nz) ! Perturbation pressure (Pascal) REAL :: qvprt (nx,ny,nz) ! Perturbation water vapor specific humidity REAL :: qc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg) REAL :: qr (nx,ny,nz) ! Rain water mixing ratio (kg/kg) REAL :: qi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg) REAL :: qs (nx,ny,nz) ! Snow mixing ratio (kg/kg) REAL :: qh (nx,ny,nz) ! Hail mixing ratio (kg/kg) REAL :: tke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2) REAL :: kmh (nx,ny,nz) ! Horizontal turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: kmv (nx,ny,nz) ! Vertical turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: ubar (nx,ny,nz) ! Base state u-velocity (m/s) REAL :: vbar (nx,ny,nz) ! Base state v-velocity (m/s) REAL :: wbar (nx,ny,nz) ! Base state w-velocity (m/s) REAL :: ptbar (nx,ny,nz) ! Base state potential temperature (K) REAL :: pbar (nx,ny,nz) ! Base state pressure (Pascal) REAL :: rhobar (nx,ny,nz) ! Base state air density (kg/m**3) REAL :: qvbar (nx,ny,nz) ! Base state water vapor specific humidity REAL :: tsoil (nx,ny,nzsoil) ! Deep soil temperature (K) REAL :: qsoil (nx,ny,nzsoil) ! Deep soil moisture REAL :: wetcanp(nx,ny) ! Canopy water amount REAL :: raing (nx,ny) ! Grid supersaturation rain REAL :: rainc (nx,ny) ! Cumulus convective rain ! !----------------------------------------------------------------------- ! ! Verification data interpolated to model grid ! !----------------------------------------------------------------------- ! REAL :: auprt (nx,ny,nz) ! Perturbation u-velocity (m/s) REAL :: avprt (nx,ny,nz) ! Perturbation v-velocity (m/s) REAL :: awprt (nx,ny,nz) ! Perturbation w-velocity (m/s) REAL :: aptprt (nx,ny,nz) ! Perturbation potential temperature (K) REAL :: apprt (nx,ny,nz) ! Perturbation pressure (Pascal) REAL :: aqvprt (nx,ny,nz) ! Perturbation water vapor specific humidity REAL :: aqc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg) REAL :: aqr (nx,ny,nz) ! Rain water mixing ratio (kg/kg) REAL :: aqi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg) REAL :: aqs (nx,ny,nz) ! Snow mixing ratio (kg/kg) REAL :: aqh (nx,ny,nz) ! Hail mixing ratio (kg/kg) REAL :: atke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2) REAL :: akmh (nx,ny,nz) ! Horizontal turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: akmv (nx,ny,nz) ! Vertical turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: aubar (nx,ny,nz) ! Base state u-velocity (m/s) REAL :: avbar (nx,ny,nz) ! Base state v-velocity (m/s) REAL :: awbar (nx,ny,nz) ! Base state w-velocity (m/s) REAL :: aptbar (nx,ny,nz) ! Base state potential temperature (K) REAL :: arhobar(nx,ny,nz) ! Base state density (kg/m**3) REAL :: apbar (nx,ny,nz) ! Base state pressure (Pascal) REAL :: aqvbar (nx,ny,nz) ! Base state water vapor specific humidity REAL :: atsoil (nx,ny,nzsoil) ! Soil temperature (K) REAL :: aqsoil (nx,ny,nzsoil) ! Soil moisture (m**3/m**3) REAL :: awetcanp(nx,ny) ! Canopy water amount REAL :: araing (nx,ny) ! Grid supersaturation rain REAL :: arainc (nx,ny) ! Cumulus convective rain ! !------------------------------------------------------------------------ ! ! Difference arrays ! !----------------------------------------------------------------------- ! REAL :: duprt (nx,ny,nz) ! perturbation x component of velocity (m/s) REAL :: dvprt (nx,ny,nz) ! perturbation y component of velocity (m/s) REAL :: dwprt (nx,ny,nz) ! perturbation vertical component of ! velocity in Cartesian coordinates (m/s) REAL :: dptprt (nx,ny,nz) ! perturbation potential temperature (K) REAL :: dpprt (nx,ny,nz) ! perturbation pressure (Pascal) REAL :: dqvprt (nx,ny,nz) ! perturbation water vapor mixing ratio (kg/kg) REAL :: dqc (nx,ny,nz) ! Cloud water mixing ratio (kg/kg) REAL :: dqr (nx,ny,nz) ! Rainwater mixing ratio (kg/kg) REAL :: dqi (nx,ny,nz) ! Cloud ice mixing ratio (kg/kg) REAL :: dqs (nx,ny,nz) ! Snow mixing ratio (kg/kg) REAL :: dqh (nx,ny,nz) ! Hail mixing ratio (kg/kg) REAL :: dtke (nx,ny,nz) ! Turbulent Kinetic Energy ((m/s)**2) REAL :: dkmh (nx,ny,nz) ! Horizontal turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: dkmv (nx,ny,nz) ! Vertical turb. mixing coef. for ! momentum. ( m**2/s ) REAL :: dtsoil (nx,ny,nzsoil) ! Soil temperature (K) REAL :: dqsoil (nx,ny,nzsoil) ! Soil moisture (m**3/m**3) REAL :: dwetcanp(nx,ny) ! Canopy water amount REAL :: draing (nx,ny) ! Grid supersaturation rain REAL :: drainc (nx,ny) ! Cumulus convective rain REAL :: tem1 (nx,ny,nz) ! A work array REAL :: tem3dsoil(nx,ny,nzsoil) ! A work array INTEGER :: ireturn, i,j,k ! !----------------------------------------------------------------------- ! ! Misc. local variables ! !----------------------------------------------------------------------- ! INTEGER :: is,js,ks,ls,ie,je,ke,le ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! is=1 js=1 ks=1 ls=1 ie=nx-1 je=ny-1 ke=nz-1 le=nzsoil-1 !----------------------------------------------------------------------- ! ! Scalars ! !----------------------------------------------------------------------- tem1=0.0 tem3dsoil = 0.0 IF(myproc == 0) PRINT *, ' ptprt: ' CALL subtrprt(nx,ny,nz, ptprt,ptbar,aptprt,aptbar,dptprt,nscale, & is,js,ks,ie,je,ke) IF(myproc == 0) PRINT *, ' pprt: ' CALL subtrprt(nx,ny,nz, pprt, pbar, apprt, apbar, dpprt,nscale, & is,js,ks,ie,je,ke) IF(myproc == 0) PRINT *, ' qvprt: ' CALL subtrprt(nx,ny,nz, qvprt,qvbar,aqvprt,aqvbar,dqvprt,nscale, & is,js,ks,ie,je,ke) ! PRINT *, ' qc: ' ! CALL subtrprt(nx,ny,nz, qc, tem1, aqc, tem1, dqc,nscale, & ! is,js,ks,ie,je,ke) ! PRINT *, ' qr: ' ! CALL subtrprt(nx,ny,nz, qr, tem1, aqr, tem1, dqr,nscale, & ! is,js,ks,ie,je,ke) ! PRINT *, ' qi: ' ! CALL subtrprt(nx,ny,nz, qi, tem1, aqi, tem1, dqi,nscale, & ! is,js,ks,ie,je,ke) ! PRINT *, ' qs: ' ! CALL subtrprt(nx,ny,nz, qs, tem1, aqs, tem1, dqs,nscale, & ! is,js,ks,ie,je,ke) ! PRINT *, ' qh: ' ! CALL subtrprt(nx,ny,nz, qh, tem1, aqh, tem1, dqh,nscale, & ! is,js,ks,ie,je,ke) ! PRINT *, ' tke: ' ! CALL subtrprt(nx,ny,nz, tke, tem1, atke, tem1, dtke,nscale, & ! is,js,ks,ie,je,ke) ! PRINT *, ' kmh: ' ! CALL subtrprt(nx,ny,nz, kmh, tem1, akmh, tem1, dkmh,nscale, & ! is,js,ks,ie,je,ke) ! PRINT *, ' kmv: ' ! CALL subtrprt(nx,ny,nz, kmv, tem1, akmv, tem1, dkmv,nscale, & ! is,js,ks,ie,je,ke) dqc=qc dqr=qr dqi=qi dqs=qs dqh=qh dtke=tke dkmh=kmh dkmv=kmv !----------------------------------------------------------------------- ! ! u wind components ! !----------------------------------------------------------------------- ie=nx IF(myproc == 0) PRINT *, ' uprt: ' CALL subtrprt(nx,ny,nz,uprt,ubar,auprt,aubar,duprt,nscale, & is,js,ks,ie,je,ke) !----------------------------------------------------------------------- ! ! v wind components ! !----------------------------------------------------------------------- ie=nx-1 je=ny IF(myproc == 0) PRINT *, ' vprt: ' CALL subtrprt(nx,ny,nz,vprt,vbar,avprt,avbar,dvprt,nscale, & is,js,ks,ie,je,ke) !----------------------------------------------------------------------- ! ! w wind components ! !----------------------------------------------------------------------- je=ny-1 ke=nz IF(myproc == 0) PRINT *, ' wprt: ' CALL subtrprt(nx,ny,nz,wprt,tem1,awprt,tem1,dwprt,nscale, & is,js,ks,ie,je,ke) !----------------------------------------------------------------------- ! ! 2-d/3-d surface (soil) variables ! !----------------------------------------------------------------------- ie=nx-1 je=ny-1 le=nzsoil ks=1 ke=1 le=1 ! ! PRINT *, ' tsoil:' ! CALL subtrprt(nx,ny,nzsoil,tsoil,tem3dsoil,atsoil,tem3dsoil,dtsoil, & ! nscale,is,js,ls,ie,je,le) ! PRINT *, ' qsoil:' ! CALL subtrprt(nx,ny,nzsoil,qsoil,tem3dsoil,aqsoil,tem3dsoil,dqsoil, & ! nscale,is,js,ls,ie,je,le) ! PRINT *, ' wetcanp:' ! CALL subtrprt(nx,ny,1,wetcanp,tem1,awetcanp,tem1,dwetcanp,nscale, & ! is,js,ks,ie,je,ke) ! PRINT *, ' raing:' ! CALL subtrprt(nx,ny,1, raing,tem1, araing,tem1, draing,nscale, & ! is,js,ks,ie,je,ke) ! PRINT *, ' rainc:' ! CALL subtrprt(nx,ny,1, rainc,tem1, arainc,tem1, drainc,nscale, & ! is,js,ks,ie,je,ke) dtsoil=tsoil dqsoil=qsoil dwetcanp=wetcanp draing=0.0 drainc=0.0 ! RETURN END SUBROUTINE prtfield ! !################################################################## !################################################################## !###### ###### !###### SUBROUTINE SUBTRPRT ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE subtrprt(nx,ny,nz, a,abar,b,bbar,c,nscale, & 6 istr,jstr,kstr,iend,jend,kend) ! !----------------------------------------------------------------------- ! ! PURPOSE: ! ! Subtracts 2 three-dimensional arrays, represented by ! means plus perturbations. ! ! AUTHOR: Keith Brewster OU School of Meteorology. Feb 1992 ! ! MODIFICATION HISTORY: ! 11 Aug 1992 (KB) changed from arps2.5 to arps3.0 ! ! !----------------------------------------------------------------------- ! ! INPUT: ! a perturbation data array ! abar mean data array ! b perturbation data array to subtract from a ! bbar mean data array to subtract from a ! ! OUTPUT: ! c difference array a-b ! (may share storage in calling program with array a or b) ! !----------------------------------------------------------------------- ! ! Variable Declarations: ! !----------------------------------------------------------------------- ! IMPLICIT NONE INCLUDE 'mp.inc' INTEGER :: nx,ny,nz ! 3 dimensions of array REAL :: nscale REAL :: a (nx,ny,nz) ! data array REAL :: abar(nx,ny,nz) ! base state of data array a REAL :: b (nx,ny,nz) ! data array to subtract from a REAL :: bbar(nx,ny,nz) ! base state of data arrya b REAL :: c (nx,ny,nz) ! difference array a-b INTEGER :: istr,jstr,kstr INTEGER :: iend,jend,kend INTEGER :: i,j,k,imid,jmid,kmid ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! ! Beginning of executable code... ! !@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ! imid=nint(0.5*FLOAT(istr+iend)) jmid=nint(0.5*FLOAT(jstr+jend)) kmid=nint(0.5*FLOAT(kstr+kend)) IF (nz < 10) kmid=kstr IF(myproc == 0) PRINT *, 'imid,jmid,kmid: ',imid,jmid,kmid ! !----------------------------------------------------------------------- ! ! Tell us about a sample input point ! !----------------------------------------------------------------------- ! IF (abs(nscale) < 1e-5) THEN IF(myproc == 0) & PRINT *, ' sample, a= ',(a(imid,jmid,kmid)+abar(imid,jmid,kmid)), & ' b= ',(b(imid,jmid,kmid)+bbar(imid,jmid,kmid)) ELSE IF(myproc == 0) & PRINT *, ' sample, a= ',(a(imid,jmid,kmid)+abar(imid,jmid,kmid)), & ' b= ',b(imid,jmid,kmid)*abs(nscale) END IF ! !----------------------------------------------------------------------- ! ! Subtraction ! !----------------------------------------------------------------------- ! DO k=kstr,kend DO j=jstr,jend DO i=istr,iend IF (abs(nscale) < 1e-5) THEN c(i,j,k)=a(i,j,k)+abar(i,j,k)-(b(i,j,k)+bbar(i,j,k)) ELSE c(i,j,k)=a(i,j,k)+b(i,j,k)*nscale END IF END DO END DO END DO ! !----------------------------------------------------------------------- ! ! Tell us about a sample output point ! !----------------------------------------------------------------------- ! IF (abs(nscale) < 1e-5) THEN IF(myproc == 0) & PRINT *, ' c= ',c(imid,jmid,kmid) ELSE IF(myproc == 0) & PRINT *, ' c= ',c(imid,jmid,kmid)+abar(imid,jmid,kmid) END IF RETURN END SUBROUTINE subtrprt !################################################################## !################################################################## !###### ###### !###### SUBROUTINE DHSLINI ###### !###### ###### !###### Developed by ###### !###### Center for Analysis and Prediction of Storms ###### !###### University of Oklahoma ###### !###### ###### !################################################################## !################################################################## ! SUBROUTINE dhslini(nx,ny,nzsoil,nstyps, & 1 tsoil,qsoil,wetcanp) IMPLICIT NONE INTEGER :: nx,ny,nzsoil,nstyps INTEGER :: i,j,k,l REAL :: tsoil (nx,ny,nzsoil,0:nstyps) ! Soil temperature (K) REAL :: qsoil (nx,ny,nzsoil,0:nstyps) ! Soil moisture (m**3/m**3) REAL :: wetcanp(nx,ny,0:nstyps) ! Canopy water amount DO k=1,nstyps DO j=1,ny DO i=1,nx DO l=1,nzsoil tsoil(i,j,l,k)=tsoil(i,j,l,0) qsoil(i,j,l,k)=qsoil(i,j,l,0) END DO wetcanp(i,j,k)=wetcanp(i,j,0) END DO END DO END DO RETURN END SUBROUTINE dhslini SUBROUTINE normalrand(iseed,n,rndn) 6,3 IMPLICIT NONE INCLUDE 'mp.inc' INTEGER :: iseed ! random seed INTEGER :: n ! dimension size INTEGER :: i, n0 REAL :: rndn (n) ! return standard normal distri. REAL, ALLOCATABLE :: work1(:),work2(:) REAL, PARAMETER :: pi=3.141592653589 REAL :: ave,var ALLOCATE (work1(n),work2(n)) n0 = n do i=1,n iseed = mod(iseed*7141+54773,259200) work1(i)=float(iseed)/259199. enddo do i=1,n iseed = mod(iseed*7141+54773,259200) work2(i)=float(iseed)/259199. enddo rndn = sqrt(-2.0*log(work1+tiny(1.0)))*cos(2.0*pi*work2) ! N(0,1) standard normal distribution ! statistics ! ave = sum(rndn)/float(n) ave = sum(rndn) IF (mp_opt > 0) THEN CALL mpsumr(ave, 1) CALL mptotali(n0) END IF ave = ave/float(n0) work1 = rndn-ave ! var = dot_product(work1,work1)/float(n-1) var = dot_product(work1,work1) IF (mp_opt > 0) CALL mpsumr(var, 1) var = var/float(n0-1) IF(myproc == 0) print *,' AVE,VAR,SQRT(AVR):',ave,var,sqrt(var) DEALLOCATE (work1,work2) RETURN END SUBROUTINE normalrand