Multi-dimentional and multivariate Time-Series forecast (RNN/LSTM) Keras
The input shape for an LSTM must be (num_samples, num_time_steps, num_features). In your example case, combining both cities as input, num_features will be 2x3=6.

If you lump all your 365 time steps into one sample, then the first dimension will be 1 - one single sample! You can also do sanity check by using the total number of data points. You have 2 cities, each with 365 time-steps and 3 features: 2x365x3= 2190 . This is obviously the same as 1x365x6 (as I said above) - so it would be a possibility (Keras will run) - but it obviously won't learn to generalise at all, only giving it one sample.

Have a look at this relevant question, which I recently answered. There I speak a little about using a rolling window (check the comments of the answer for more info). That will buy you more samples if you need them.

If you want to train a single model with data for both cities as input, then making predictions for both cities at each time-step is as simple as defining a final Dense layer, which outputs 2 units. Your validation/test data must then of course contain a tuple of (city1, city2).

A perhaps more sophisticated way to approach this would be to create datasets on a single-city basis, then train several sub-models on each city individually (say for 5 layers), then Merge/Concatenate them and put several further layers on top. This will mean you are combining the learnt features of each city, which are in turn being combined to a higher level of abstraction. Here is the first image I got from a search engine, which sketches the idea.