I have demonstrated in my own tests that explicit division of work among threads, so as to use static scheduling, performs much better than dynamic on multi-CPU platforms. On a single CPU with shared last level cache, one of schedule(dynamic,2) or (guided,2) or the like should work. If you must choose among static, dynamic, and guided, static may work OK in spite of some thread taking double the average work; evidently then, dynamic or guided pinned to a single CPU socket may also work.
Thread private arrays significantly increase various overheads involved in starting or completing parallel regions. Even with thread private scalars, there appears to be a relatively low number of distinct private variables which can work effectively with Intel OpenMP. The compiler's -parallel option may be superior to OpenMP when an excessive number of privates is involved.
I've also made some comparisons with Cilk+. The latter becomes more competitive against OpenMP as dynamic scheduling comes into play. Cilk+ institutes a great divide between simple cases which are more straightforward than OpenMP and more complex reductions where the opposite is true.

Petros,

I suggest you create a module that contains

A(:,:,:),IPVT(:,:),Ty(:,:),Tx(:,:),Bx(:,:),By(:,:),TB(:,:),x(:,:),f(:,:), Cy, Cx

And then USE that module within your two subroutines

!$omp parallel do reduction(+:counter)
doi=1,num_domains
callcalculate_domain(i)
callsolve_domain(i)
counter=counter+1
enddo
!$omp end parallel do

Jim Dempsey

>>Is there any benefit instead of explicitly declaring them as shared?

You won't be passing all of the arguments on each call to the two inner subroutines.
All that is required to be passed is the index i.

You did not indicate if the called routine takes a assumed shape arrays. My guess is yes, because of the way the arguments are declared on the call. You are making 10,000 calls. You did not list the extent of the other dimensions on the called arrays. If they are relatively small, then the overhead to construct and pass the array descriptors may be a significant portion of the loop. If you have VTune available, then you can measure the time spent making the calls (as opposed to the time spent in what was called).

Jim Dempsey

Your first post had num_domains at 10000 (not 1000), and num_elem at 40 therefore

doi=1,num_domains
callcalculate_domain(A(:,:,i),IPVT(:,i),Ty(:,i),Tx(:,i),Bx(:,i),By(:,i),TB(:,i),x(:,i),f(:,i),Cy,Cx)
callsolve_domain(A(:,:,i),IPVT(:,i),Ty(:,i),Tx(:,i),Bx(:,i),By(:,i),TB(:,i),x(:,i),f(:,i))
counter=counter+1
enddo

would:

loop:
reserve stack space for array descriptors (hardly any time)
construct array descriptors A(:,:,i),IPVT(:,i),Ty(:,i),Tx(:,i),Bx(:,i),By(:,i),TB(:,i),x(:,i),f(:,i),Cy(:),Cx(:) each possibly 20 (machine) words to initialize (x11 array args) IOW at least 220 machine words to write
call calculate_domain (hardly any time
iterate 40 times on all but A
return
reserve stack space for array descriptors (hardly any time)
construct array descriptors
A(:,:,i),IPVT(:,i),Ty(:,i),Tx(:,i),Bx(:,i),By(:,i),TB(:,i),x(:,i),f(:,i)
each possibly 20 (machine) words to initialize (x9 array args) IOW at
least 180 machine words to write
call calculate_domain (hardly any time
iterate 40 times on all but A
return
end loop

Depending on what is done in each subroutine, with 40 iterations (assumption) the time to construct the array descriptors could exceed the processing time within the called subroutine.

If these subroutines can be inlined, then the compiler optimization might possibly remove some of the redundancy. e.g. construct the array descriptors once and/or pushing (some of) the array descriptor construction out of the loop. I cannot assess as if the optimization will do this. However, you can easily do this by using the module technique.

Jim Dempsey