Implementing the projected atomic orbital (PAO) localization in Gaussian

I am currently attempting to learn how to setup the active space for performing CASSCF calculations. 


The textbook I am using as a reference is:


*Roos, Björn O. "Multiconfigurational quantum chemistry." Theory and Applications of Computational Chemistry. 2005. 725-764.*


My specific problem is discussed in Section 13.1.1. I am working with Gaussian 2016.


I would like to generate a checkpoint file with orbitals localized as atomic orbitals in order to use them as a starting point when selecting my active space. 
What I mean is the following...


I am working with BODIPY, with 11 active orbitals (the PI framework).


From a normal SCF I obtain the typical MOs, these would be three of my active orbitals:


[![enter image description here](https://i.stack.imgur.com/UqiZd.jpg)](https://i.stack.imgur.com/UqiZd.jpg)


If I include the option Pop=SaveNBOs, I obtain localized Natural Bonding Orbitals. These would be three of my NBO active orbitals:


[![enter image description here](https://i.stack.imgur.com/8GANm.jpg)](https://i.stack.imgur.com/8GANm.jpg)


I have found that I can use IOp 4/9 to control the localization method to some extent (for example, "IOp(4/9=20212)" implements the Pipek-Mezey localization). However, after searching for a while, I still don't see how to generate the localized atomic orbitals. Could any of you suggest how to achieve this? I would then be working with the 11 pz orbitals as my active space. 


Thanks!


