search

geoschem / geos-chem

GEOS-Chem "Science Codebase" repository. Contains GEOS-Chem science routines, run directory generation scripts, and interface code. This repository is used as a submodule within the GCClassic and GCHP wrappers, as well as in other modeling contexts (external ESMs).
http://geos-chem.org
Other
168 stars 165 forks source link

[BUG/ISSUE] Potential error in toms_mod.F90 in GEOS-Chem 13 #591

Closed yantosca closed 3 years ago

yantosca commented 3 years ago

Tess Carter wrote:

I am attaching a terminal readout of the issue, which seems to stem from a State_Chm initialization or lack thereof of TOMS1_O3_COL and TOMS2_O3_COL.

Based on my understanding of the default options in the model, we should be using online met fields for ozone columns and not even getting here in the code, so I have 2 questions.

  1. By commenting out that part of the if statement (attached toms_mod.F90 file), the model compiles fine. Is there a scientific problem with doing that?
  2. Is there an easy fix to this issue? We were surprised to see it come up. We haven’t modified anything else in the code.

Thanks in advance, Tess

Screen Shot 2021-01-12 at 9 29 48 PM
yantosca commented 3 years ago

Thanks for bringing this issue to our attention.

From the error output you sent, I believe that State_Chm should be declared with INTENT(INOUT) instead of INTENT(IN).

Ideally the model should be picking up the TO3 from the met fields, unless there is something going on that I’m not aware of. I am going to look into this now and hope to have some clarity on this soon.

yantosca commented 3 years ago

What I tried

So I did a test with gfortran 10.2.0. I made sure that State_Chm was INTENT(IN) in COMPUTE_OVERHEAD_O3. Then I put some debug print in both COMPUTE_OVERHEAD_O3 (in toms_mod.F90), and again in GET_OVERHEAD_O3_FOR_FASTJ (in Interfaces/GCClassic/main.F90). I ran it expecting an error (since we are trying to modify State_Chm, which is read-only. But instead, I got this:

==============================================================================
---> DATE: 2019/07/01  UTC: 00:10  X-HRS:      0.166667
---> DATE: 2019/07/01  UTC: 00:20  X-HRS:      0.333333
### in compute_overhead_o3: use_o3_from_met
### sum to3_daily, sum_to3:    979499.37290954590        979499.37290954590     
 ### SUM( State_Chm%TO3_DAILY ) in main:    979499.37290954590     
     - DO_STRAT_CHEM: Linearized strat chemistry at 2019/07/01 00:20
     - LINOZ_CHEM3: Doing LINOZ
---> DATE: 2019/07/01  UTC: 00:30  X-HRS:      0.500000
---> DATE: 2019/07/01  UTC: 00:40  X-HRS:      0.666667
### in compute_overhead_o3: use_o3_from_met
### sum to3_daily, sum_to3:    979499.37290954590        979499.37290954590     
 ### SUM( State_Chm%TO3_DAILY ) in main:    979499.37290954590     
...

so that meant that State_Chm%TO3_DAILY was getting set to the proper value of column overhead O3. Huh?

Explanation

So this is my shot at an explanation. The routine COMPUTE_OVERHEAD_O3 is called from COMPUTE_OVERHEAD_O3_FOR_FASTJ, which is an internal routine in Interfaces/GCClassic/main.F90. The calling sequence is this:

       !=====================================================================
       !                  ***** C H E M I S T R Y *****
       !=====================================================================
       IF ( notDryRun ) THEN
          IF ( Input_Opt%useTimers ) THEN
             CALL Timer_Start( "All chemistry", RC )
          ENDIF

          ! Get the overhead column O3 for use with FAST-J
          ! NOTE: Move to CHEMISTRY section.  This now has to come after
          ! the call to HEMCO emissions driver EMISSIONS_RUN. (bmy, 3/20/15)
          CALL Get_Overhead_O3_For_FastJ( Input_Opt, State_Grid, State_Met )

and the subroutine is this:

!------------------------------------------------------------------------------
!                  GEOS-Chem Global Chemical Transport Model                  !
!------------------------------------------------------------------------------
!BOP
!
! !IROUTINE: get_overhead_o3_for_fastj
!
! !DESCRIPTION: Internal subroutine GET\_OVERHEAD\_O3\_FOR\_FASTJ
!\\
!\\
! !INTERFACE:
!
  SUBROUTINE Get_Overhead_O3_For_FastJ( Input_Opt, State_Grid, State_Met )
!
! !USES:
!
    USE Input_Opt_Mod,      ONLY : OptInput
    USE State_Grid_Mod,     ONLY : GrdState
    USE State_Met_Mod,      ONLY : MetState
!
! !INPUT PARAMETERS:
!
    TYPE(OptInput), INTENT(IN)    :: Input_Opt
    TYPE(GrdState), INTENT(IN)    :: State_Grid
    TYPE(MetState), INTENT(IN)    :: State_Met
!
! !REMARKS:
!  This routine makes use of variables declared in above in the main program
!  (which are visible in all sub-programs below the CONTAINS statement).
!                                                                             .
!  The original code was done in FAST-J routine "set_prof.F", but has been
!  split off to facilitate development of the grid-independent model.
!
! !REVISION HISTORY:
!  07 Mar 2012 - R. Yantosca - Initial version
!  See https://github.com/geoschem/geos-chem for complete history
!EOP
!------------------------------------------------------------------------------
!BOC
!
! !LOCAL VARIABLES:
!
    ! FAST-J is only used for fullchem and offline aerosol, skip otherwise
    IF ( ITS_A_FULLCHEM_SIM .or. ITS_AN_AEROSOL_SIM  ) THEN

       ! Only execute this if we are doing chemistry
       ! and if it we are at a chemistry timestep
       IF ( LCHEM .and. ITS_TIME_FOR_CHEM() ) THEN

          ! Get the overhead O3 column for FAST-J.  Take either the
          ! TOMS O3 data or the column O3 directly from the met fields
          CALL Compute_Overhead_O3( Input_Opt, State_Grid, State_Chm, DAY, &
                                    Input_Opt%USE_O3_FROM_MET,  &
                                    State_Met%TO3, RC )
          print*, '### SUM( State_Chm%TO3_DAILY ) in main: ', &
               sum(state_chm%to3_daily)
       ENDIF
    ENDIF

  END SUBROUTINE Get_Overhead_O3_For_FastJ

Note that we don't pass State_Chm via an argument-to-dummy-variable (as we do for Input_Opt, State_Met, State_Grid). In that case, GET_OVERHEAD_O3_FOR_FASTJ is using the global variable State_Chm, which is declared at the top of the main.F90 program. I believe that the compiler that I used then passes State_Chm to GET_OVERHEAD_O3_FOR_FASTJ as a read/write variable (i.e. INTENT(INOUT)). This might be because gfortran 10.2 skews more to the the F2018 standards than to strict F90. But older compilers (like your ifort12) might still throw an error.

Solution

Therefore, I believe the most robust thing to do is to pass State_Chm as INTENT(INOUT) to GET_OVERHEAD_O3_FOR_FASTJ:

  SUBROUTINE Get_Overhead_O3_For_FastJ( Input_Opt,  State_Chm,               &
                                        State_Grid, State_Met, RC            )
!
! !USES:
!
    USE Input_Opt_Mod,      ONLY : OptInput
    USE State_Chm_Mod,      ONLY : ChmState
    USE State_Grid_Mod,     ONLY : GrdState
    USE State_Met_Mod,      ONLY : MetState
!
! !INPUT PARAMETERS:
!
    TYPE(OptInput), INTENT(IN)    :: Input_Opt
    TYPE(GrdState), INTENT(IN)    :: State_Grid
    TYPE(MetState), INTENT(IN)    :: State_Met
!
! !INPUT/OUTPUT PARAMETERS:
!
    TYPE(ChmState), INTENT(INOUT) :: State_Chm
!
! OUTPUT PARAMETERS:
!
    INTEGER,        INTENT(OUT)   :: RC
   .. etc...

and then also declare State_Chm as INTENT(INOUT) in GET_OVERHEAD_O3 in toms_mod.F90:

  SUBROUTINE COMPUTE_OVERHEAD_O3( Input_Opt, State_Grid,      State_Chm,     &
                                  DAY,       USE_O3_FROM_MET, TO3,           &
                                  RC                                        )
!
! !USES:
!
    USE HCO_Calc_Mod,     ONLY : Hco_EvalFld
    USE HCO_State_GC_Mod, ONLY : HcoState
    USE Input_Opt_Mod,    ONLY : OptInput
    USE State_Grid_Mod,   ONLY : GrdState
    USE State_Chm_Mod,    ONLY : ChmState
!
! !INPUT PARAMETERS:
!
    TYPE(OptInput), INTENT(IN)    :: Input_Opt             ! Input Options
    TYPE(GrdState), INTENT(IN)    :: State_Grid            ! Grid State
    INTEGER,        INTENT(IN)    :: DAY                   ! Day of month
    LOGICAL,        INTENT(IN)    :: USE_O3_FROM_MET       ! Use TO3 rom met?
    REAL(fp),       INTENT(IN)    :: TO3(State_Grid%NX, &
                                         State_Grid%NY)    ! Met TO3 [Dobsons]
!
! !INPUT/OUTPUT PARAMETERS
!
    TYPE(ChmState), INTENT(INOUT) :: State_Chm             ! Chemistry State
!
! !OUTPUT PARAMETERS: 
!
    INTEGER,        INTENT(OUT)   :: RC                    ! Success/failure?
!
... etc. ...

We can also pass back the RC flag so that we can stop the code in the main routine in case of error. That should be agnostic of compiler version.

I will push this to the dev/13.0.0-final branch, so that it goes into the final 13.0.0 release.

Good catch!!

yantosca commented 3 years ago

This issue has been resolved with commit 1630a809ea822f39318b0e2b35badb2e71c38697. It will ship with the 13.0.0 final release version.