First comment: you should be using real data, rather than simulated data, as we are not very good at modelling genetic variation and recombination at a genome-wide scale. A way that this can be done with GCTA is mentioned in the PLINK manual.
I'll do some calculations to estimate the size of data you need to store:
Assume that this plink-formatted file is in tped format (one row per SNP), and the entire genome is diploid. That means that a single SNP/individual pair will take up three bytes. You want 30,000 cases and 70,000 controls, so each line from each sample will take up 300kB (plus a bit of change for the header). With 1 thousand SNPs, that's 300MB per sample; with 1 million SNPs, that's 300GB per sample; with 10 million SNPs, that's 3TB per sample; with 40 million SNPs, that's 12 TB for one sample.
Reading the plink manual, I see that the BED format is a two-bit compression for each SNP/individual pair, so one byte per 4 individuals, or about 12x compression (i.e. 1TB for one sample).
Okay, that's a lot. It's well within a 256TB storage limit, but might hit file size limits depending on the file system.
If there is a data generation issue with PLINK, it might be worth trying (as you suggest) splitting the input for simulation then merging afterwards. My guess is that doing the simulation on fewer SNPs (rather than fewer individuals) would work better. Plink has a
--bmerge option for merging the compressed binary file formats.
I usually approach problems like these by working my way up from something that I know works. Try first merging two datasets of 1000 SNPs. If that works, try larger datasets. Rinse and repeat until the computer crashes, you get bored, or you reach your goal.