diff --git a/README.standalone b/README.standalone new file mode 100644 index 00000000..15e06515 --- /dev/null +++ b/README.standalone @@ -0,0 +1,33 @@ +Instructions on how to build and run the stand alone ca_global code on hera +1. Copy enitre directory strucure + cd /scratch2/BMC/gsienkf/Philip.Pegion + cp -r NEMSfv3gfs_develop + +2. Optional Compile full model: (I alredy built is, so you should not have to) + cd /tests/ + ./compile.sh $PWD/../FV3 hera.intel '32BIT=Y HYDRO=N' 1 NO NO + +4. compile standalone model + cd ../stochastic_physics + # my loaded modules: + 1) intel/18.0.5.274 3) netcdf/4.7.0 5) contrib 7) anaconda/latest + 2) impi/2018.0.4 4) hdf5/1.10.5 6) grads/2.2.1 + + + compile_standalone_ca OR compile_standalone + +5. There is an input.nml which is currently set to run at C96, and the INPUT directory points to + C96 grid files. If you want to run at C384, you will need to rm INPUT and link in the C384 + directory (ln -sf INPUT_C384 INPUT) and make npx and npy 385 in input.nml + +5. get an interactive session + + salloc --x11=first -q debug -t 0:30:00 --nodes=1 -A -n 6 + +6. run job + srun standalone_ca OR standalone_stochy + +7. post-process job to 360x180 grid (script is on hera in ~Philip.Pegion) + run_fregrid.sh (OR run_fregrid_c384.sh) + + output file is ca_out.nc diff --git a/compile_standalone b/compile_standalone index f947baaa..11daf3e0 100755 --- a/compile_standalone +++ b/compile_standalone @@ -38,6 +38,4 @@ if [ $compile_all -eq 1 ];then $FC ${FLAGS64} stochastic_physics.F90 ar rv libstochastic_physics.a *.o fi -#$FC ${FLAGS64} get_stochy_pattern.F90 -#ar rv libstochastic_physics.a *.o $FC -traceback -real-size 64 -qopenmp -o standalone_stochy standalone_stochy.F90 ${INCS} -I/apps/netcdf/4.7.0/intel/18.0.5.274/include -L. -lstochastic_physics -L../FV3/atmos_cubed_sphere -lfv3core -L../FMS/FMS_INSTALL -lfms -L../FV3/gfsphysics -lgfsphys -L/scratch2/NCEPDEV/nwprod/NCEPLIBS/compilers/intel/18.0.5.274/lib -lsp_v2.0.3_d -L/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -Wl,-rpath,/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -lesmf -L/apps/netcdf/4.7.0/intel/18.0.5.274/lib -lnetcdff -lnetcdf diff --git a/compile_standalone_ca b/compile_standalone_ca new file mode 100755 index 00000000..cce7e8ce --- /dev/null +++ b/compile_standalone_ca @@ -0,0 +1,17 @@ +#!/bin/sh +compile_all=1 +if [ $compile_all -eq 1 ];then + rm -f *.i90 *.i *.o *.mod lib*a + mpif90 -C -traceback -real-size 64 -c ../FV3/gfsphysics/physics/machine.F + mpif90 -C -traceback -real-size 64 -c ../FV3/gfsphysics/physics/mersenne_twister.f + mpif90 -C -traceback -I../FMS/include -c ../FMS/platform/platform.F90 + mpif90 -C -traceback -nowarn -DGFS_PHYS -I../FMS/include -c ../FMS/constants/constants.F90 + mpif90 -C -traceback -I../FV3/atmos_cubed_sphere -I../FMS/include -I../FMS/fv3gfs -c fv_control_stub.F90 + mpif90 -C -traceback -I../FV3/atmos_cubed_sphere -I../FMS/include -I../FMS/fv3gfs -c atmosphere_stub.F90 + mpif90 -C -traceback -c -real-size 64 standalone_stochy_module.F90 + mpif90 -C -traceback -I. -real-size 64 -c plumes.F90 + mpif90 -DSTOCHY_UNIT_TEST -real-size 64 -C -traceback -I../FMS/fv3gfs -I../FMS/FMS_INSTALL -I../FV3/atmos_cubed_sphere -c update_ca.F90 + mpif90 -DSTOCHY_UNIT_TEST -real-size 64 -C -traceback -I../FMS/fv3gfs -I../FMS/FMS_INSTALL -I../FV3/atmos_cubed_sphere -c cellular_automata_global.F90 + ar rv libcellular_automata.a *.o +fi +mpif90 -traceback -real-size 64 -qopenmp -o standalone_ca standalone_ca.F90 -I../FV3/atmos_cubed_sphere -I../FMS/FMS_INSTALL -I/apps/netcdf/4.7.0/intel/18.0.5.274/include -L. -lcellular_automata -L../FV3/atmos_cubed_sphere -lfv3core -L../FMS/FMS_INSTALL -lfms -L../FV3/gfsphysics -lgfsphys -L/scratch2/NCEPDEV/nwprod/NCEPLIBS/compilers/intel/18.0.5.274/lib -lsp_v2.0.3_d -L/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -Wl,-rpath,/scratch1/NCEPDEV/nems/emc.nemspara/soft/esmf/8.0.0bs48-intel18.0.5.274-impi2018.0.4-netcdf4.6.1/lib -lesmf -L/apps/netcdf/4.7.0/intel/18.0.5.274/lib -lnetcdff -lnetcdf diff --git a/compns_stochy.F90 b/compns_stochy.F90 index d116cf23..a87b32ab 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -44,8 +44,10 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) character(len=*), intent(in) :: input_nml_file(sz_nml) character(len=64), intent(in) :: fn_nml real, intent(in) :: deltim - real tol + real tol,l_min + real :: rerth,circ,tmp_lat integer k,ios + integer,parameter :: four=4 ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! @@ -58,6 +60,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) namelist /nam_sfcperts/nsfcpert,pertz0,pertshc,pertzt,pertlai, & ! mg, sfcperts pertvegf,pertalb,iseed_sfc,sfc_tau,sfc_lscale,sppt_land + rerth =6.3712e+6 ! radius of earth (m) tol=0.01 ! tolerance for calculations ! spectral resolution defintion ntrunc=-999 @@ -120,7 +123,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) ! length scale. skeb_varspect_opt = 0 sppt_logit = .false. ! logit transform for sppt to bounded interval [-1,+1] - fhstoch = -999.0 ! forecast hour to dump random patterns + fhstoch = -999.0 ! forecast interval (in hours) to dump random patterns stochini = .false. ! true= read in pattern, false=initialize from seed #ifdef INTERNAL_FILE_NML @@ -207,6 +210,33 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret) pertlai(1) > 0 .OR. pertvegf(1) > 0 .OR. pertalb(1) > 0) THEN do_sfcperts=.true. ENDIF +!calculate ntrunc if not supplied + if (ntrunc .LT. 1) then + if (me==0) print*,'ntrunc not supplied, calculating' + circ=2*3.1415928*rerth ! start with lengthscale that is circumference of the earth + l_min=circ + do k=1,5 + if (sppt(k).GT.0) l_min=min(sppt_lscale(k),l_min) + if (shum(k).GT.0) l_min=min(shum_lscale(k),l_min) + if (skeb(k).GT.0) l_min=min(skeb_lscale(k),l_min) + enddo + !ntrunc=1.5*circ/l_min + ntrunc=circ/l_min + if (me==0) print*,'ntrunc calculated from l_min',l_min,ntrunc + endif + ! ensure lat_s is a mutiple of 4 with a reminader of two + ntrunc=INT((ntrunc+1)/four)*four+2 + if (me==0) print*,'NOTE ntrunc adjusted for even nlats',ntrunc + +! set up gaussian grid for ntrunc if not already defined. + if (lon_s.LT.1 .OR. lat_s.LT.1) then + lat_s=ntrunc*1.5+1 + lon_s=lat_s*2+4 +! Grid needs to be larger since interpolation is bi-linear + lat_s=lat_s*2 + lon_s=lon_s*2 + if (me==0) print*,'gaussian grid not set, defining here',lon_s,lat_s + endif ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ! All checks are successful. diff --git a/get_stochy_pattern.F90 b/get_stochy_pattern.F90 index c157317d..f682ee09 100644 --- a/get_stochy_pattern.F90 +++ b/get_stochy_pattern.F90 @@ -81,6 +81,7 @@ subroutine get_random_pattern_fv3(rpattern,npatterns,& enddo call mp_reduce_sum(workg,lonf,latg) + ! interpolate to cube grid allocate(rslmsk(lonf,latg)) diff --git a/standalone_ca.F90 b/standalone_ca.F90 new file mode 100644 index 00000000..8c84b7b4 --- /dev/null +++ b/standalone_ca.F90 @@ -0,0 +1,265 @@ +program standalone_stochy_new + +use standalone_stochy_module + +use fv_mp_mod, only: mp_start, domain_decomp +use atmosphere_stub_mod, only: Atm,atmosphere_init_stub +use fv_arrays_mod, only: fv_atmos_type +use fv_control_mod, only: setup_pointers +!use mpp_domains +use mpp_mod, only: mpp_set_current_pelist,mpp_get_current_pelist,mpp_init,mpp_pe,mpp_npes ,mpp_declare_pelist +use mpp_domains_mod, only: mpp_broadcast_domain,MPP_DOMAIN_TIME,mpp_domains_init ,mpp_domains_set_stack_size +use fms_mod, only: fms_init +!use time_manager_mod, only: time_type +use xgrid_mod, only: grid_box_type +use netcdf + + +implicit none +type(GFS_control_type) :: Model +integer, parameter :: nlevs=64 +integer :: ntasks,fid,ct +integer :: nthreads,omp_get_num_threads +integer :: ncid,xt_dim_id,yt_dim_id,time_dim_id,xt_var_id,yt_var_id,time_var_id,ca_out_id +integer :: ca1_id,ca2_id,ca3_id +character*1 :: strid +type(GFS_grid_type),allocatable :: Grid(:) +type(GFS_diag_type),allocatable :: Diag(:) +type(GFS_statein_type),allocatable :: Statein(:) +type(GFS_coupling_type),allocatable :: Coupling(:) +include 'mpif.h' +include 'netcdf.inc' +real(kind=4) :: ts,undef + +integer :: cres,blksz,nblks,ierr,my_id,i,j,nx2,ny2,nx,ny,id +integer,target :: npx,npy +integer :: ng,layout(2),io_layout(2),commID,grid_type,ntiles +integer :: halo_update_type = 1 +logical,target :: nested +integer :: pe,npes,stackmax=4000000 + +real(kind=4),allocatable,dimension(:,:) :: workg +real(kind=4),allocatable,dimension(:) :: grid_xt,grid_yt +real(kind=8),pointer ,dimension(:,:) :: area +type(grid_box_type) :: grid_box +!type(time_type) :: Time ! current time +!type(time_type) :: Time_step ! atmospheric time step. +!type(time_type) :: Time_init ! reference time. +!---cellular automata control parameters +integer :: nca !< number of independent cellular automata +integer :: nlives !< cellular automata lifetime +integer :: ncells !< cellular automata finer grid +real :: nfracseed !< cellular automata seed probability +integer :: nseed !< cellular automata seed frequency +logical :: do_ca !< cellular automata main switch +logical :: ca_sgs !< switch for sgs ca +logical :: ca_global !< switch for global ca +logical :: ca_smooth !< switch for gaussian spatial filter +logical :: isppt_deep !< switch for combination with isppt_deep. OBS! Switches off SPPT on other tendencies! +logical :: isppt_pbl +logical :: isppt_shal +logical :: pert_flux +logical :: pert_trigger +integer :: iseed_ca !< seed for random number generation in ca scheme +integer :: nspinup !< number of iterations to spin up the ca +integer :: ca_amplitude +real :: nthresh !< threshold used for perturbed vertical velocity + +NAMELIST /gfs_physics_nml/ nca, ncells, nlives, nfracseed,nseed, nthresh, & + do_ca,ca_sgs, ca_global,iseed_ca,ca_smooth,isppt_pbl,isppt_shal,isppt_deep,nspinup,& + pert_trigger,pert_flux,ca_amplitude + +! default values +nca = 1 +ncells = 5 +nlives = 10 +nfracseed = 0.5 +nseed = 100000 +iseed_ca = 0 +nspinup = 1 +do_ca = .false. +ca_sgs = .false. +ca_global = .false. +ca_smooth = .false. +isppt_deep = .false. +isppt_shal = .false. +isppt_pbl = .false. +pert_trigger = .false. +pert_flux = .false. +ca_amplitude = 500. +nthresh = 0.0 + +! open namelist file +open (unit=565, file='input.nml', READONLY, status='OLD', iostat=ierr) +read(565,gfs_physics_nml) +close(565) +! define stuff +ng=3 ! ghost region +undef=9.99e+20 + +! initialize fms +call fms_init() +call mpp_init() +call fms_init +my_id=mpp_pe() + +call atmosphere_init_stub (grid_box, area) +!define domain +isd=Atm(1)%bd%isd +ied=Atm(1)%bd%ied +jsd=Atm(1)%bd%jsd +jed=Atm(1)%bd%jed +isc=Atm(1)%bd%isc +iec=Atm(1)%bd%iec +jsc=Atm(1)%bd%jsc +jec=Atm(1)%bd%jec +nx=Atm(1)%npx-1 +ny=Atm(1)%npy-1 +allocate(workg(nx,ny)) + +! for this simple test, nblocks = ny, blksz=ny +nblks=ny +blksz=nx +nthreads = omp_get_num_threads() +Model%me=my_id + +Model%nca = nca +Model%ncells = ncells +Model%nlives = nlives +Model%nfracseed = nfracseed +Model%nseed = nseed +Model%iseed_ca = iseed_ca +Model%nspinup = nspinup +Model%do_ca = do_ca +Model%ca_sgs = ca_sgs +Model%ca_global = ca_global +Model%ca_smooth = ca_smooth +Model%isppt_deep = isppt_deep +Model%isppt_pbl = isppt_pbl +Model%isppt_shal = isppt_shal +Model%pert_flux = pert_flux +Model%pert_trigger = pert_trigger +Model%nthresh = nthresh + +! setup GFS_init parameters + +!define model grid + +allocate(grid_xt(nx),grid_yt(ny)) +do i=1,nx + grid_xt(i)=i +enddo +do i=1,ny + grid_yt(i)=i +enddo + +!setup GFS_coupling +allocate(Diag(nblks)) +allocate(Coupling(nblks)) +allocate(Statein(nblks)) +write(strid,'(I1.1)') my_id+1 +fid=30+my_id +ierr=nf90_create('ca_out.tile'//strid//'.nc',cmode=NF90_CLOBBER,ncid=ncid) +ierr=NF90_DEF_DIM(ncid,"grid_xt",nx,xt_dim_id) +ierr=NF90_DEF_DIM(ncid,"grid_yt",ny,yt_dim_id) +ierr=NF90_DEF_DIM(ncid,"time",NF90_UNLIMITED,time_dim_id) + !> - Define the dimension variables. +ierr=NF90_DEF_VAR(ncid,"grid_xt",NF90_FLOAT,(/ xt_dim_id /), xt_var_id) +ierr=NF90_PUT_ATT(ncid,xt_var_id,"long_name","T-cell longitude") +ierr=NF90_PUT_ATT(ncid,xt_var_id,"cartesian_axis","X") +ierr=NF90_PUT_ATT(ncid,xt_var_id,"units","degrees_E") +ierr=NF90_DEF_VAR(ncid,"grid_yt",NF90_FLOAT,(/ yt_dim_id /), yt_var_id) +ierr=NF90_PUT_ATT(ncid,yt_var_id,"long_name","T-cell latitude") +ierr=NF90_PUT_ATT(ncid,yt_var_id,"cartesian_axis","Y") +ierr=NF90_PUT_ATT(ncid,yt_var_id,"units","degrees_N") +ierr=NF90_DEF_VAR(ncid,"time",NF90_FLOAT,(/ time_dim_id /), time_var_id) +ierr=NF90_PUT_ATT(ncid,time_var_id,"long_name","time") +ierr=NF90_PUT_ATT(ncid,time_var_id,"units","hours since 2014-08-01 00:00:00") +ierr=NF90_PUT_ATT(ncid,time_var_id,"cartesian_axis","T") +ierr=NF90_PUT_ATT(ncid,time_var_id,"calendar_type","JULIAN") +ierr=NF90_PUT_ATT(ncid,time_var_id,"calendar","JULIAN") +!ierr=NF90_DEF_VAR(ncid,"ca_out",NF90_FLOAT,(/xt_dim_id, yt_dim_id ,time_dim_id/), ca_out_id) +!ierr=NF90_PUT_ATT(ncid,ca_out_id,"long_name","random pattern") +!ierr=NF90_PUT_ATT(ncid,ca_out_id,"units","None") +!ierr=NF90_PUT_ATT(ncid,ca_out_id,"missing_value",undef) +!ierr=NF90_PUT_ATT(ncid,ca_out_id,"_FillValue",undef) +!ierr=NF90_PUT_ATT(ncid,ca_out_id,"cell_methods","time: point") +ierr=NF90_DEF_VAR(ncid,"ca1",NF90_FLOAT,(/xt_dim_id, yt_dim_id ,time_dim_id/), ca1_id) +ierr=NF90_PUT_ATT(ncid,ca1_id,"long_name","random pattern") +ierr=NF90_PUT_ATT(ncid,ca1_id,"units","None") +ierr=NF90_PUT_ATT(ncid,ca1_id,"missing_value",undef) +ierr=NF90_PUT_ATT(ncid,ca1_id,"_FillValue",undef) +ierr=NF90_PUT_ATT(ncid,ca1_id,"cell_methods","time: point") +ierr=NF90_DEF_VAR(ncid,"ca2",NF90_FLOAT,(/xt_dim_id, yt_dim_id ,time_dim_id/), ca2_id) +ierr=NF90_PUT_ATT(ncid,ca2_id,"long_name","random pattern") +ierr=NF90_PUT_ATT(ncid,ca2_id,"units","None") +ierr=NF90_PUT_ATT(ncid,ca2_id,"missing_value",undef) +ierr=NF90_PUT_ATT(ncid,ca2_id,"_FillValue",undef) +ierr=NF90_PUT_ATT(ncid,ca2_id,"cell_methods","time: point") +ierr=NF90_DEF_VAR(ncid,"ca3",NF90_FLOAT,(/xt_dim_id, yt_dim_id ,time_dim_id/), ca3_id) +ierr=NF90_PUT_ATT(ncid,ca3_id,"long_name","random pattern") +ierr=NF90_PUT_ATT(ncid,ca3_id,"units","None") +ierr=NF90_PUT_ATT(ncid,ca3_id,"missing_value",undef) +ierr=NF90_PUT_ATT(ncid,ca3_id,"_FillValue",undef) +ierr=NF90_PUT_ATT(ncid,ca3_id,"cell_methods","time: point") +ierr=NF90_ENDDEF(ncid) +ierr=NF90_PUT_VAR(ncid,xt_var_id,grid_xt) +ierr=NF90_PUT_VAR(ncid,yt_var_id,grid_yt) +! allocate diagnostics +DO i =1,nblks + allocate(Diag(i)%ca_out(blksz)) + allocate(Diag(i)%ca_deep(blksz)) + allocate(Diag(i)%ca_turb(blksz)) + allocate(Diag(i)%ca_shal(blksz)) + allocate(Diag(i)%ca_rad(blksz)) + allocate(Diag(i)%ca_micro(blksz)) + allocate(Diag(i)%ca1(blksz)) + allocate(Diag(i)%ca2(blksz)) + allocate(Diag(i)%ca3(blksz)) +! allocate coupling + allocate(Coupling(i)%cape(blksz)) + allocate(Coupling(i)%ca_out(blksz)) + allocate(Coupling(i)%ca_deep(blksz)) + allocate(Coupling(i)%ca_turb(blksz)) + allocate(Coupling(i)%ca_shal(blksz)) + allocate(Coupling(i)%ca_rad(blksz)) + allocate(Coupling(i)%ca_micro(blksz)) + allocate(Coupling(i)%ca1(blksz)) + allocate(Coupling(i)%ca2(blksz)) + allocate(Coupling(i)%ca3(blksz)) +! allocate coupling + allocate(Statein(i)%pgr(blksz)) + allocate(Statein(i)%qgrs(blksz,nlevs,1)) + allocate(Statein(i)%vvl(blksz,nlevs)) + allocate(Statein(i)%prsl(blksz,nlevs)) +ENDDO +ct=1 +do i=1,600 + ts=i/8.0 ! hard coded to write out hourly based on a 450 second time-step + call cellular_automata_global(i-1, Statein, Coupling, Diag, & + nblks, Model%levs, Model%nca, Model%ncells, & + Model%nlives, Model%nfracseed, Model%nseed, & + Model%nthresh, Model%ca_global, Model%ca_sgs, & + Model%iseed_ca, Model%ca_smooth, Model%nspinup, & + blksz) + if (mod(i,8).EQ.0) then + do j=1,ny + workg(:,j)=Diag(j)%ca1(:) + enddo + ierr=NF90_PUT_VAR(ncid,ca1_id,workg,(/1,1,ct/)) + do j=1,ny + workg(:,j)=Diag(j)%ca2(:) + enddo + ierr=NF90_PUT_VAR(ncid,ca2_id,workg,(/1,1,ct/)) + do j=1,ny + workg(:,j)=Diag(j)%ca3(:) + enddo + ierr=NF90_PUT_VAR(ncid,ca3_id,workg,(/1,1,ct/)) + ierr=NF90_PUT_VAR(ncid,time_var_id,ts,(/ct/)) + ct=ct+1 + if (my_id.EQ.0) write(6,fmt='(a,i5,4f6.3)') 'ca=',i,Diag(1)%ca1(1:4) + endif +enddo +!close(fid) +ierr=NF90_CLOSE(ncid) +end diff --git a/standalone_stochy.F90 b/standalone_stochy.F90 index a9fcbcbb..75e9ce74 100644 --- a/standalone_stochy.F90 +++ b/standalone_stochy.F90 @@ -70,7 +70,10 @@ program standalone_stochy integer :: halo_update_type = 1 real :: dx,dy,pi,rd,cp logical,target :: nested -integer :: nlunit,pe,npes,stackmax=4000000 +logical :: write_this_tile +integer :: nargs,ntile_out,nlunit,pe,npes,stackmax=4000000 +character*80 :: fname +character*1 :: ntile_out_str real(kind=4),allocatable,dimension(:,:) :: workg real(kind=4),allocatable,dimension(:,:,:) :: workg3d @@ -86,7 +89,13 @@ program standalone_stochy skeb,skeb_tau,skeb_vdof,skeb_lscale,iseed_skeb,skeb_vfilt,skeb_diss_smooth, & skeb_sigtop1,skeb_sigtop2,skebnorm,sppt_sigtop1,sppt_sigtop2,& shum_sigefold,spptint,shumint,skebint,skeb_npass,use_zmtnblck,new_lscale - +write_this_tile=.false. +ntile_out_str='0' +nargs=iargc() +if (nargs.EQ.1) then + call getarg(1,ntile_out_str) +endif +read(ntile_out_str,'(I1.1)') ntile_out open (unit=nlunit, file='input.nml', READONLY, status='OLD') read(nlunit,nam_stochy) close(nlunit) @@ -167,7 +176,7 @@ program standalone_stochy allocate(Grid(i)%xlat(blksz)) allocate(Grid(i)%xlon(blksz)) enddo -do j=1,ny +do j=1,nblks Grid(j)%xlat(:)=Init_parm%xlat(:,j) Grid(j)%xlon(:)=Init_parm%xlon(:,j) enddo @@ -181,9 +190,15 @@ program standalone_stochy !setup GFS_coupling allocate(Coupling(nblks)) call init_stochastic_physics(Model, Init_parm, ntasks, nthreads) +call get_outfile(fname) write(strid,'(I1.1)') my_id+1 +if (ntile_out.EQ.0) write_this_tile=.true. +if ((my_id+1).EQ.ntile_out) write_this_tile=.true. +print*,trim(fname)//'.tile'//strid//'.nc',write_this_tile +if (write_this_tile) then fid=30+my_id -ierr=nf90_create('workg.tile'//strid//'.nc',cmode=NF90_CLOBBER,ncid=ncid) +!ierr=nf90_create(trim(fname)//'.tile'//strid//'.nc',cmode=NF90_CLOBBER,ncid=ncid) +ierr=nf90_create(trim(fname)//'.tile'//strid//'.nc',cmode=NF90_CLOBBER,ncid=ncid) ierr=NF90_DEF_DIM(ncid,"grid_xt",nx,xt_dim_id) ierr=NF90_DEF_DIM(ncid,"grid_yt",ny,yt_dim_id) if (Model%do_skeb)ierr=NF90_DEF_DIM(ncid,"p_ref",nlevs,zt_dim_id) @@ -246,6 +261,7 @@ program standalone_stochy if (Model%do_skeb)then ierr=NF90_PUT_VAR(ncid,zt_var_id,pressl) endif +endif do i=1,nblks if (Model%do_sppt)allocate(Coupling(i)%sppt_wts(blksz,nlevs)) @@ -253,11 +269,12 @@ program standalone_stochy if (Model%do_skeb)allocate(Coupling(i)%skebu_wts(blksz,nlevs)) if (Model%do_skeb)allocate(Coupling(i)%skebv_wts(blksz,nlevs)) enddo -do i=1,40 +do i=1,200 Model%kdt=i ts=i/4.0 call run_stochastic_physics(Model, Grid, Coupling, nthreads) if (Model%me.EQ.0) print*,'sppt_wts=',i,Coupling(1)%sppt_wts(1,20) + if (write_this_tile) then if (Model%do_sppt)then do j=1,ny workg(:,j)=Coupling(j)%sppt_wts(:,20) @@ -285,6 +302,17 @@ program standalone_stochy ierr=NF90_PUT_VAR(ncid,varid4,workg3d,(/1,1,1,i/)) endif ierr=NF90_PUT_VAR(ncid,time_var_id,ts,(/i/)) + endif enddo -ierr=NF90_CLOSE(ncid) +if (write_this_tile) ierr=NF90_CLOSE(ncid) +end +subroutine get_outfile(fname) +use stochy_namelist_def +character*80,intent(out) :: fname +character*4 :: s_ntrunc,s_lat,s_lon + write(s_ntrunc,'(I4)') ntrunc + write(s_lat,'(I4)') lat_s + write(s_lon,'(I4)') lon_s + fname=trim('workg_T'//trim(adjustl(s_ntrunc))//'_'//trim(adjustl(s_lon))//'x'//trim(adjustl(s_lat))) + return end diff --git a/stochastic_physics.F90 b/stochastic_physics.F90 index 33f849b1..9b899020 100644 --- a/stochastic_physics.F90 +++ b/stochastic_physics.F90 @@ -1,13 +1,20 @@ +!>@brief The module 'stochastic_physics' is for initialization and running of +!! the stochastic physics random pattern generators module stochastic_physics implicit none private -public :: init_stochastic_physics, run_stochastic_physics +public :: init_stochastic_physics +public :: run_stochastic_physics contains +!>@brief The subroutine 'init_stochastic_physics' initializes the stochastic +!!pattern genertors +!>@details It reads the stochastic physics namelist (nam_stoch and nam_sfcperts) +!allocates and polulates the necessary arrays subroutine init_stochastic_physics(Model, Init_parm, ntasks, nthreads) use fv_mp_mod, only : is_master use stochy_internal_state_mod @@ -17,7 +24,7 @@ subroutine init_stochastic_physics(Model, Init_parm, ntasks, nthreads) use stochy_gg_def,only : colrad_a use stochy_namelist_def use physcons, only: con_pi -use spectral_layout_mod,only:me,ompthreads +use spectral_layout_mod,only:me,ompthreads,nodes use mpp_mod #ifdef STOCHY_UNIT_TEST use standalone_stochy_module, only: GFS_control_type, GFS_init_type @@ -35,7 +42,7 @@ subroutine init_stochastic_physics(Model, Init_parm, ntasks, nthreads) integer :: iret real*8 :: PRSI(Model%levs),PRSL(Model%levs),dx real, allocatable :: skeb_vloc(:) -integer :: k,kflip,latghf,nodes,blk,k2 +integer :: k,kflip,latghf,blk,k2 character*2::proc ! Set/update shared variables in spectral_layout_mod @@ -191,7 +198,8 @@ subroutine init_stochastic_physics(Model, Init_parm, ntasks, nthreads) !print *,'done with init_stochastic_physics' end subroutine init_stochastic_physics - +!>@brief The subroutine 'run_stochastic_physics' updates the random patterns if +!!necessary subroutine run_stochastic_physics(Model, Grid, Coupling, nthreads) use fv_mp_mod, only : is_master use stochy_internal_state_mod @@ -229,7 +237,7 @@ subroutine run_stochastic_physics(Model, Grid, Coupling, nthreads) maxlen = maxval(Model%blksz(:)) ! check to see if it is time to write out random patterns -if (Model%phour .EQ. fhstoch) then +if (fhstoch.GE. 0 .AND. MOD(Model%phour,fhstoch) .EQ. 0) then write(STRFH,FMT='(I6.6)') nint(Model%phour) sfile='stoch_out.F'//trim(STRFH) call dump_patterns(sfile) diff --git a/stochy_data_mod.F90 b/stochy_data_mod.F90 index c1065ba5..6291f8dc 100644 --- a/stochy_data_mod.F90 +++ b/stochy_data_mod.F90 @@ -59,15 +59,15 @@ subroutine init_stochdata(nlevs,delt,input_nml_file,fn_nml,nlunit,iret) levs=nlevs iret=0 - if(is_master()) print*,'in init stochdata' call compns_stochy (me,size(input_nml_file,1),input_nml_file(:),fn_nml,nlunit,delt,iret) + if(is_master()) print*,'in init stochdata',nodes,lat_s if ( (.NOT. do_sppt) .AND. (.NOT. do_shum) .AND. (.NOT. do_skeb) .AND. (.NOT. do_sfcperts) ) return - if (nodes.GE.lat_s/2) then - lat_s=(int(nodes/12)+1)*24 - lon_s=lat_s*2 - ntrunc=lat_s-2 - if (is_master()) print*,'WARNING: spectral resolution is too low for number of mpi_tasks, resetting lon_s,lat_s,and ntrunc to',lon_s,lat_s,ntrunc - endif +! if (nodes.GE.lat_s/2) then +! lat_s=(int(nodes/12)+1)*24 +! lon_s=lat_s*2 +! ntrunc=lat_s-2 +! if (is_master()) print*,'WARNING: spectral resolution is too low for number of mpi_tasks, resetting lon_s,lat_s,and ntrunc to',lon_s,lat_s,ntrunc +! endif call initialize_spectral(gis_stochy, iret) if (iret/=0) return allocate(noise_e(len_trie_ls,2),noise_o(len_trio_ls,2)) diff --git a/sumfln_stochy.f b/sumfln_stochy.f index 91a751fa..cbbfd23f 100644 --- a/sumfln_stochy.f +++ b/sumfln_stochy.f @@ -18,7 +18,8 @@ subroutine sumfln_stochy(flnev,flnod,lat1s,plnev,plnod, use machine use spectral_layout_mod, only : num_parthds_stochy => ompthreads !or : use fv_mp_mod ? - use mpp_mod, only: mpp_npes, mpp_alltoall, mpp_get_current_pelist + use mpp_mod, only: mpp_pe,mpp_npes, mpp_alltoall, + & mpp_get_current_pelist implicit none ! @@ -70,11 +71,11 @@ subroutine sumfln_stochy(flnev,flnod,lat1s,plnev,plnod, ! real(kind=4),dimension(2*nvars*ls_dim*workdim*nodes):: ! & work1dr,work1ds real(kind=4),pointer:: work1dr(:),work1ds(:) - integer, dimension(jcap+1) :: kpts, kptr, sendcounts, recvcounts, + integer, dimension(nodes) :: kpts, kptr, sendcounts, recvcounts, & sdispls ! integer ierr,ilat,ipt_ls, lmax,lval,i,jj,lonl,nv - integer node,nvar,arrsz + integer node,nvar,arrsz,my_pe integer ilat_list(nodes) ! for OMP buffer copy ! ! statement functions @@ -93,6 +94,7 @@ subroutine sumfln_stochy(flnev,flnod,lat1s,plnev,plnod, num_threads = min(num_parthds_stochy,nvars) nvar_thread_max = (nvars+num_threads-1)/num_threads npes = mpp_npes() + my_pe=mpp_pe() allocate(pelist(0:npes-1)) call mpp_get_current_pelist(pelist) kpts = 0 @@ -187,7 +189,6 @@ subroutine sumfln_stochy(flnev,flnod,lat1s,plnev,plnod, do node = 1, nodes - 1 ilat_list(node+1) = ilat_list(node) + lats_nodes(node) end do - !$omp parallel do private(node,jj,ilat,lat,ipt_ls,nvar,kn,n2) do node=1,nodes do jj=1,lats_nodes(node)