running GCHP adjoint, code crashes

[kaysuselj]

Your name

kay suselj

Your affiliation

JPL

What happened? What did you expect to happen?

I am running the adjoint version of the GCHP using CO2 w/ CMS-Flux emissions setup.
The forward model run is executed successfully with checkpoint files being written out.

I have problems with adjoint:

In src/GCHP_GridComp/GEOSChem_GridComp/geos-chem/Interfaces/GCHP/Chem_GridCompMod.F90 subroutine Initialize_ lines 1650+:
DO I = 1, nFlds
N = State_Chm%Map_Advect(I)
ThisSpc => State_Chm%SpcData(N)%Info
Int2Adj(I)%Name = TRIM(ThisSpc%Name)

seem to create array that contains adjoint variables, for all species.

The execution stalls at I=2 with forrtl: severe (408): fort: (2): Subscript #1 of the array MAP_ADVECT has value 3 which is greater than the upper bound of 2
In my case nFlds=2 and Map_Advect(2)=3.

Any suggestion on whether this is a bug in the code or with my setup would be appreciated.

Thanks.

What are the steps to reproduce the bug?

run adjoint version of GCHP.

Please attach any relevant configuration and log files.

No response

What GCHP version were you using?

14.5.1

What environment were you running GCHP on?

Other (please explain below)

What compiler and version were you using?

pleiades

What MPI library and version were you using?

mpi

Will you be addressing this bug yourself?

Yes, but I will need some help

Additional information

No response

[lizziel]

Hi @kaysuselj, I will take a look at the code. Could you also post your GEOS-Chem log file? If you could run with verbose true in geoschem_config.yml that would be great.

[kaysuselj]

Thanks @lizziel . Attached is the gchp.log. A couple of comments:

• lines that start with Kay ... - this are lines I added to the code to print a couple of additional outputs, you can ignore them.
• I see that I get the unit conversion error earlier in the code, but then the code seems to keep running. Is that the problem?

gchp.log

[lizziel]

Thanks @kaysuselj. I think you should comment out the lines in input_mod.F90 that are causing the early errors. It is trying to print values from geoschem_config.yml to the file. We have not run into an issue with those prints before. What compiler are you using? Regardless, best remove that problem since it is the not the main one.

I do not see what the Map_Advect issue is from just the log. In your prints could you write the names of the variables you are printing? That might help me determine what I am seeing. I think your debug code is going in the right direction.

[kaysuselj]

Thanks, Lizzie.

Hm, I commented a few lines in the input_mod.F90 that were causing the first error:
WRITE( 6,110 ) 'GAMMA HO2 : ', Input_Opt%GAMMA_HO2
WRITE( 6,120 ) 'Absolute AR threshold : ', Input_Opt%AUTOREDUCE_THRESHOL
WRITE( 6, 120 ) 'Ice supersaturation req. : ', & ! ((Input_Opt%P_ICE_SUPERSAT-1)*1.e+2_fp)

My understanding is that these were only the diagnostic outputs, and did not modify the model run.

I recompiled and rerun the gchp with these changes. To my surprise, the code runs past the previous point of failure.
However, now the problem seems to be with reading the checkpoint files (assigning the variables to the CONC pointer in Chem_GridCompMod.F90) on line: State_Chm%Species(IND)%Conc(:,:,:) = & Ptr3D_R8(:,:,State_Grid%NZ:1:-1)

Commenting the lines above, have I changed something else in the code?

---

1. Compilers: using ifort on pleiades, the attached file gchp.sh are the setting that I use.
2. I am attaching my version of the Chem_GridComp.F90, it is probably easiest for you to see what my write-outs are.
3. Attached is also the last version of gchp.log

gchp.sh.txt

Chem_GridCompMod.F90.txt

gchp.log

---

Thanks a lot, Kay

[lizziel]

Hi @kaysuselj, I'm glad that worked. I suspect that the issue with the print was corrupting a core which resulted in being messed up later on.

The new issue is because the adjoint code has not been tested since an update to point the State_Chm species arrays directly at the internal state. Copying is no longer done. If you look a little further down in the code from where your code fails you will see where the concentration is pointed to the MAPL values. It is here:
https://github.com/geoschem/geos-chem/blob/86cd586260e4f2a9db7efd06cea19a5fe0bdc9e8/Interfaces/GCHP/Chem_GridCompMod.F90#L2597-L2603

I also noticed that the section right after that (here) shows that the adjoint species concentration arrays still do copies and the old way of indexing. I recommend updating that for consistency with the rest of the model. Let me know if you need any tips for that. It basically would mimic what we do for State_Chm%Species%Conc except for adjoint species arrays.

[kaysuselj]

Hi @lizziel thanks, I just want to make sure i completely understand what the code does before making changes. From your comments, I understand there are 2 issues:

The code crash in this block. My understanding here is that:

1. line 2560 we update the internal state from the checkpoint file
2. line 2569: Ptr3D_R8 points to the concentration of species that we are working on. Is that necessary to do? My understanding is that the Int2Spc already points to the new internal state (i.e. the one modified from the checkpoint file).
3. If 2 is true, that means on line 2573, I need to replace Ptr3D_R8(:,:,State_Grid%NZ:1:-1) with Int2Spc(IND)%Internal(:,:,State_Grid%NZ:1:-1)?

Concentrations of adjoint variables -- i see that in the old version the adjoint concentrations are 4D arrays with the last dimension representing the species, and we should change this. I will probably need to modify this in all routines where State_Chm%SpeciesAdj is called. I would like to address the first crash before doing that.

[lizziel]

Hi @kaysuselj, the reason the code is crashing is the concentration array is a pointer that is not yet pointing at anything. To put the internal state value into it you need to use the code at https://github.com/geoschem/geos-chem/blob/86cd586260e4f2a9db7efd06cea19a5fe0bdc9e8/Interfaces/GCHP/Chem_GridCompMod.F90#L2597-L2603.

My understanding here is that:

1. line 2560 we update the internal state from the checkpoint file

I am not familiar with that adjoint subroutine. For regular GCHP we do not have to call anything special to get the internal state from the checkpoint since it is automatic in MAPL. I looked over the subroutine quickly but it has not documentation. Could you add some at some point?

2. line 2569: Ptr3D_R8 points to the concentration of species that we are working on. Is that necessary to do? My understanding is that the Int2Spc already points to the new internal state (i.e. the one modified from the checkpoint file).
   Correct, you should use Int2Spc, and you need to point the Conc arrays at it rather than use equals. Note this is done later, so you need to do some updating of the ADJOINT blocks to avoid doing it twice.

3. If 2 is true, that means on line 2573, I need to replace Ptr3D_R8(:,:,State_Grid%NZ:1:-1) with Int2Spc(IND)%Internal(:,:,State_Grid%NZ:1:-1)?
   Yes, but again, you need to update the CCP blocks so as not to do this twice.

[kaysuselj]

Hi @lizziel sorry for delay.

I got the adjoint working :)

Here are a couple of observations:

1. The issue with the crashes related to writeouts from input_mod.F90 (running GCHP adjoint, code crashes #483 (comment)) seem to be present only if I compile the code with the debug option (flags -DCMAKE_BUILD_TYPE=Debug).
2. The current version only works for the CO2 w/ CMS-Flux emissions option. There are still some issues when I try to run the full carbon simulation with the 4 species.

Here are the required modifications to make adjoint run

• GCHP code: https://github.com/kaysuselj/GCHP
  
• geos_chem https://github.com/kaysuselj/geos-chem_adjoint
• HEMCO https://github.com/kaysuselj/HEMCO_adjoint

Is this the best way for me to suggest changes and if so, can you check them and include in the new version of repos?

Thanks.

[lizziel]

Hi Kay, could you make a draft pull request for each of the repositories with changes? You can do this by clicking "Create pull request" on your GitHub forks. This will make it easier to review your changes and bring them in. Let me know if you need any help with this.

[kaysuselj]

Done. Please let me know if you have any questions. Kay