Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add fixes for TOMAS simulations in GCHP #2136

Merged
merged 6 commits into from
Feb 6, 2024
Merged

Add fixes for TOMAS simulations in GCHP #2136

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
c3c6b71
Bug fixes for TOMAS
BettyCroft Jan 12, 2024
dea155d
revert to original createRunDir - only need to remove -P for WashU ma…
BettyCroft Jan 12, 2024
1fd5c26
Bug fix: Remove hardcoded TOMAS species from GCHP's HISTORY.rc.fullchem
msulprizio Jan 17, 2024
2d0cd6b
Merge pull request #2116 for branch for branch 'tomas/dev' into dev/1…
msulprizio Jan 24, 2024
f4ed069
Update restart file path for GCHP TOMAS simulations
msulprizio Jan 30, 2024
ca65457
Fix precision in code modified by TOMAS updates
msulprizio Feb 5, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
- Increse requested time limits in GCHP integration tests (compile 2h30m, run 5h)
- Changed CO2 concentration used in RRTMG to be modifiable in geoschem_config.yml
- Changed water vapor used in RRTMG to match to tracer field at all altitudes
- Updated restart file path for GCHP TOMAS simulations

### Fixed
- Fixed bug in stratospheric aerosols optical depths passed to Fast-JX
Expand Down
2 changes: 1 addition & 1 deletion GeosCore/carbon_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4968,7 +4968,7 @@ SUBROUTINE EMISSCARBONTOMAS( Input_Opt, State_Chm, State_Grid, State_Met, &
REAL*4 :: TEMPTMS ! calculated from State_Met
REAL*4 :: PRES ! calculated from State_Met
REAL(fp) :: OC2OM = 1.8d0
LOGICAL :: SGCOAG = .True.
LOGICAL :: SGCOAG = .FALSE. ! bc,jrp turn off subgrid coag 18/12/23
INTEGER :: L, K, EMTYPE
INTEGER :: ii=53, jj=29
INTEGER :: origUnit
Expand Down
8 changes: 6 additions & 2 deletions GeosCore/isorropiaII_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -572,14 +572,18 @@ SUBROUTINE DO_ISORROPIAII( Input_Opt, State_Chm, State_Diag, &
!-----------------------------------------------

IF (N == 1) THEN
IF (Spc(id_SALAAL)%Conc(I,J,L) .GT. CONMIN) THEN
! IF (Spc(id_SALAAL)%Conc(I,J,L) .GT. CONMIN) THEN !bc 08/01/24 for floating point exception
IF (Spc(id_SALAAL)%Conc(I,J,L) .GT. CONMIN .and. &
Spc(id_SALA)%Conc(I,J,L) .GT. CONMIN) THEN
AlkR = Spc(id_SALAAL)%Conc(I,J,L) / Spc(id_SALA)%Conc(I,J,L)
AlkR = MAX( (1.0_fp-AlkR), CONMIN)
ELSE
AlkR = 1.0_fp
ENDIF
ELSE
IF (Spc(id_SALCAL)%Conc(I,J,L) .GT. CONMIN) THEN
! IF (Spc(id_SALCAL)%Conc(I,J,L) .GT. CONMIN) THEN !bc 08/01/24 for floating point exception
IF (Spc(id_SALCAL)%Conc(I,J,L) .GT. CONMIN .and. &
Spc(id_SALC)%Conc(I,J,L) .GT. CONMIN) THEN
AlkR = Spc(id_SALCAL)%Conc(I,J,L) / Spc(id_SALC)%Conc(I,J,L)
AlkR = MAX( (1.0_fp-AlkR), CONMIN)
ELSE
Expand Down
23 changes: 12 additions & 11 deletions GeosCore/sulfate_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1000,7 +1000,7 @@ SUBROUTINE SRCSF30( Input_Opt, State_Grid, State_Met, State_Chm, TC2, RC )

INTEGER :: N_TRACERS
INTEGER :: IBINS
LOGICAL :: ERRORSWITCH, SGCOAG = .TRUE.
LOGICAL :: ERRORSWITCH, SGCOAG = .FALSE. ! bc,jrp - turn off subgridcoag 18/12/23
INTEGER :: FLAG, ERR
logical :: pdbug !(temporary) win, 10/24/07
!integer :: ii, jj, ll
Expand Down Expand Up @@ -1045,17 +1045,18 @@ SUBROUTINE SRCSF30( Input_Opt, State_Grid, State_Met, State_Chm, TC2, RC )
AREA = HcoState%Grid%AREA_M2%Val(:,:)
!AREA3D(:,:,1) = AREA(:,:)
!AREA3D(:,:,2) = AREA(:,:)


! comment out if shut off subgridcoag bc,14/12/23
! Define subgrid coagulation timescale (win, 10/28/08)
IF ( TRIM(State_Grid%GridRes) == '4.0x5.0' ) THEN
TSCALE = 10.*3600. ! 10 hours
ELSE IF ( TRIM(State_Grid%GridRes) == '2.0x2.5' ) THEN
TSCALE = 5.*3600.
ELSE IF ( TRIM(State_Grid%GridRes) == '0.5x0.625' ) THEN
TSCALE = 1.*3600.
ELSE IF ( TRIM(State_Grid%GridRes) == '0.25x0.3125' ) THEN
TSCALE = 0.5*3600.
ENDIF
!IF ( TRIM(State_Grid%GridRes) == '4.0x5.0' ) THEN
! TSCALE = 10.0_fp*3600.0_fp ! 10 hours
!ELSE IF ( TRIM(State_Grid%GridRes) == '2.0x2.5' ) THEN
TSCALE = 5.0_fp*3600.0_fp
!!ELSE IF ( TRIM(State_Grid%GridRes) == '0.5x0.625' ) THEN
! TSCALE = 1.0_fp*3600.0_fp
!ELSE IF ( TRIM(State_Grid%GridRes) == '0.25x0.3125' ) THEN
! TSCALE = 0.5_fp*3600.0_fp
!ENDIF
msulprizio marked this conversation as resolved.
Show resolved Hide resolved

! Point to species array
TC1 => State_Chm%Species
Expand Down
251 changes: 203 additions & 48 deletions GeosCore/tomas_mod.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -524,41 +524,41 @@ SUBROUTINE AEROPHYS( Input_Opt, State_Chm, State_Grid, State_Met, &
printneg = .FALSE.

! determine ion rate
!ionrate = 10. ! set as constant now !!jrp

IF ( TRIM(State_Grid%GridRes) == '4.0x5.0' ) THEN

if( (pres / 100.) .lt. cplevels(9) ) then
ionrate = cosmic_ions(i,j,9) * boxmass / boxvol
if(State_Met%FRCLND(I,J) .gt. 0.2) then
ionrate = ionrate + soil_ions(9)
endif
elseif((pres/100.) .gt. cplevels(1)) then
ionrate = cosmic_ions(i,j,1) * boxmass / boxvol
if( State_Met%FRCLND(I,J) .gt. 0.2 ) then
ionrate = ionrate + soil_ions(1)
endif
else
lev=2
do while (pres / 100. .lt. cplevels(lev))
lev=lev+1
enddo
weight=( cplevels( lev - 1 ) - pres / 100. ) / &
( cplevels( lev - 1 ) - cplevels(lev) )
ionrate=( cosmic_ions(i,j,lev ) * weight + &
cosmic_ions(i,j,lev-1) * (1.e+0_fp - weight) ) &
* boxmass / boxvol
if( State_Met%FRCLND(I,J) .gt. 0.2) then
ionrate=ionrate + ( soil_ions( lev ) * weight + &
soil_ions( lev-1 ) * (1.e+0_fp-weight) )
endif
endif

ELSE
ionrate = 0.e+0_fp
ENDIF

if(ionrate .le. 1.501) ionrate = 1.501
ionrate = 10.0_fp ! set as constant now !!jrp, bc 18/12/23 and comment out below

! IF ( TRIM(State_Grid%GridRes) == '4.0x5.0' ) THEN
!
! if( (pres / 100.) .lt. cplevels(9) ) then
! ionrate = cosmic_ions(i,j,9) * boxmass / boxvol
! if(State_Met%FRCLND(I,J) .gt. 0.2) then
! ionrate = ionrate + soil_ions(9)
! endif
! elseif((pres/100.) .gt. cplevels(1)) then
! ionrate = cosmic_ions(i,j,1) * boxmass / boxvol
! if( State_Met%FRCLND(I,J) .gt. 0.2 ) then
! ionrate = ionrate + soil_ions(1)
! endif
! else
! lev=2
! do while (pres / 100. .lt. cplevels(lev))
! lev=lev+1
! enddo
! weight=( cplevels( lev - 1 ) - pres / 100. ) / &
! ( cplevels( lev - 1 ) - cplevels(lev) )
! ionrate=( cosmic_ions(i,j,lev ) * weight + &
! cosmic_ions(i,j,lev-1) * (1.e+0_fp - weight) ) &
! * boxmass / boxvol
! if( State_Met%FRCLND(I,J) .gt. 0.2) then
! ionrate=ionrate + ( soil_ions( lev ) * weight + &
! soil_ions( lev-1 ) * (1.e+0_fp-weight) )
! endif
! endif
!
! ELSE
! ionrate = 0.e+0_fp
! ENDIF
!
! if(ionrate .le. 1.501) ionrate = 1.501
msulprizio marked this conversation as resolved.
Show resolved Hide resolved

!print*,'i',i,'j',j,'l',l,'ionrate',ionrate

Expand Down Expand Up @@ -731,8 +731,16 @@ SUBROUTINE AEROPHYS( Input_Opt, State_Chm, State_Grid, State_Met, &
ERR_IND(2) = J
ERR_IND(3) = L
ERR_IND(4) = 0
! IF (SPINUP(14.0) .and. Gcout(jc) /= Gcout(jc) ) THEN
IF( SPINUP(14.0) .AND. IT_IS_NAN( Gcout(jc) ) ) THEN
Gcout(jc) = 0.0e+0_fp ! reset Nan to zero during spinup, bc 18/12/23
print*,'Reset Gcout NaN to zero at ',I,J,L
ELSEIF ( SPINUP(14.0) .AND. .not. IT_IS_FINITE( Gcout(jc) ) ) THEN
Gcout(jc) = 0.0e+0_fp ! reset Inf to zero during spinup, bc 18/12/23
print*,'Reset Gcout Inf to zero at ',I,J,L
ELSE
call check_value( Gcout(jc), ERR_IND, ERR_VAR, ERR_MSG )

ENDIF
!if( IT_IS_FINITE(Gcout(jc))) then
! print *,'xxxxxxxxx Found Inf in Gcout xxxxxxxxxxxxxx'
! print *,'Location ',I,J,L, 'comp',jc
Expand Down Expand Up @@ -3891,7 +3899,7 @@ SUBROUTINE SOACOND( MSOA, I, J, L, BOXVOL, TEMPTMS, PRES, BOXMASS,&

! For SOACOND warnings
INTEGER, SAVE :: SOACOND_WARNING_CT = -1
INTEGER, PARAMETER :: SOACOND_WARNING_MAX = 20
INTEGER, PARAMETER :: SOACOND_WARNING_MAX = 2 !bc 02/01/24 shut off too many warnings

!=================================================================
! SOACOND begins here
Expand Down Expand Up @@ -6044,6 +6052,150 @@ SUBROUTINE AERODIAG( PTYPE, I, J, L, Nk, Nkd, Mk, Mkd, BOXMASS, &
ENDDO
ENDIF

IF ( PTYPE == 2 ) THEN
DO K=1,IBINS
State_Diag%TomasCOAGmass(I,J,L,K) = 0.e+0_fp
State_Diag%TomasCOAGnumber(I,J,L,K) = 0.e+0_fp
DO JS = 1, ICOMP-IDIAG
IF ( State_Diag%Archive_TomasCOAGmass ) THEN
S = State_Diag%Map_TomasCOAGmass%id2slot(K)
IF ( S > 0 ) THEN
State_Diag%TomasCOAGmass(I,J,L,K) = &
State_Diag%TomasCOAGmass(I,J,L,K) + (MK(K,JS) - MKD(K,JS))/ &
DTCHEM / BOXMASS ! kg/kg air/sec
ENDIF
ENDIF
ENDDO
IF ( State_Diag%Archive_TomasCOAGnumber ) THEN
S = State_Diag%Map_TomasCOAGnumber%id2slot(K)
IF ( S > 0 ) THEN
State_Diag%TomasCOAGnumber(I,J,L,K) = &
State_Diag%TomasCOAGnumber(I,J,L,K) + (NK(K) - NKD(K))/ &
DTCHEM / BOXMASS ! no./kg air/sec
ENDIF
ENDIF
ENDDO
ENDIF

IF ( PTYPE == 3 ) THEN
DO K=1,IBINS
State_Diag%TomasNUCLmass(I,J,L,K) = 0.e+0_fp
State_Diag%TomasNUCLnumber(I,J,L,K) = 0.e+0_fp
DO JS = 1, ICOMP-IDIAG
IF ( State_Diag%Archive_TomasNUCLmass ) THEN
S = State_Diag%Map_TomasNUCLmass%id2slot(K)
IF ( S > 0 ) THEN
State_Diag%TomasNUCLmass(I,J,L,K) = &
State_Diag%TomasNUCLmass(I,J,L,K) + (MK(K,JS) - MKD(K,JS))/ &
DTCHEM / BOXMASS ! kg/kg air/sec
ENDIF
ENDIF
ENDDO

IF ( State_Diag%Archive_TomasNUCLnumber ) THEN
S = State_Diag%Map_TomasNUCLnumber%id2slot(K)
IF ( S > 0 ) THEN
State_Diag%TomasNUCLnumber(I,J,L,K) = &
State_Diag%TomasNUCLnumber(I,J,L,K) + (NK(K) - NKD(K))/ &
DTCHEM / BOXMASS ! no./kg air/sec
ENDIF
ENDIF
ENDDO
ENDIF

IF ( PTYPE == 4 ) THEN
DO K=1,IBINS
State_Diag%TomasAQOXmass(I,J,L,K) = 0.e+0_fp
State_Diag%TomasAQOXnumber(I,J,L,K) = 0.e+0_fp
DO JS = 1, ICOMP-IDIAG
IF ( State_Diag%Archive_TomasAQOXmass ) THEN
S = State_Diag%Map_TomasAQOXmass%id2slot(K)
IF ( S > 0 ) THEN
State_Diag%TomasAQOXmass(I,J,L,K) = &
State_Diag%TomasAQOXmass(I,J,L,K) + (MK(K,JS) - MKD(K,JS))/ &
DTCHEM / BOXMASS ! kg/kg air/sec
ENDIF
ENDIF
ENDDO
IF ( State_Diag%Archive_TomasAQOXnumber ) THEN
S = State_Diag%Map_TomasAQOXnumber%id2slot(K)
IF ( S > 0 ) THEN
State_Diag%TomasAQOXnumber(I,J,L,K) = &
State_Diag%TomasAQOXnumber(I,J,L,K) + (NK(K) - NKD(K))/ &
DTCHEM / BOXMASS ! no./kg air/sec \

ENDIF
ENDIF
ENDDO
ENDIF

IF ( PTYPE == 5 ) THEN
DO K=1,IBINS
State_Diag%TomasMNFIXmass(I,J,L,K) = 0.e+0_fp
State_Diag%TomasMNFIXnumber(I,J,L,K) = 0.e+0_fp
DO JS = 1, ICOMP-IDIAG
IF ( State_Diag%Archive_TomasMNFIXmass ) THEN
S = State_Diag%Map_TomasMNFIXmass%id2slot(K)
IF ( S > 0 ) THEN
State_Diag%TomasMNFIXmass(I,J,L,K) = &
State_Diag%TomasMNFIXmass(I,J,L,K) + (MK(K,JS) - MKD(K,JS))/ &
DTCHEM / BOXMASS ! kg/kg air/sec
ENDIF
ENDIF
ENDDO
IF ( State_Diag%Archive_TomasMNFIXnumber ) THEN
S = State_Diag%Map_TomasMNFIXnumber%id2slot(K)
IF ( S > 0 ) THEN
State_Diag%TomasMNFIXnumber(I,J,L,K) = &
State_Diag%TomasMNFIXnumber(I,J,L,K) + (NK(K) - NKD(K))/ &
DTCHEM / BOXMASS ! no./kg air/sec
ENDIF
ENDIF
ENDDO
ENDIF

IF ( PTYPE == 6 ) THEN
DO K=1,IBINS
State_Diag%TomasSOAmass(I,J,L,K) = 0.e+0_fp
State_Diag%TomasSOAnumber(I,J,L,K) = 0.e+0_fp
DO JS = 1, ICOMP-IDIAG
IF ( State_Diag%Archive_TomasSOAmass ) THEN
S = State_Diag%Map_TomasSOAmass%id2slot(K)
IF ( S > 0 ) THEN
State_Diag%TomasSOAmass(I,J,L,K) = &
State_Diag%TomasSOAmass(I,J,L,K) + (MK(K,JS) - MKD(K,JS))/ &
DTCHEM / BOXMASS ! kg/kg air/sec
ENDIF
ENDIF
ENDDO
IF ( State_Diag%Archive_TomasSOAnumber ) THEN
S = State_Diag%Map_TomasSOAnumber%id2slot(K)
IF ( S > 0 ) THEN
State_Diag%TomasSOAnumber(I,J,L,K) = &
State_Diag%TomasSOAnumber(I,J,L,K) + (NK(K) - NKD(K))/ &
DTCHEM / BOXMASS ! no./kg air/sec
ENDIF
ENDIF
ENDDO
ENDIF

IF ( PTYPE == 7 ) THEN
DO K=1,IBINS
State_Diag%TomasNUCRATEnumber(I,J,L,K) = 0.e+0_fp
IF ( State_Diag%Archive_TomasNUCRATEnumber ) THEN
S = State_Diag%Map_TomasNUCRATEnumber%id2slot(K)
IF ( S > 0 ) THEN
State_Diag%TomasNUCRATEnumber(I,J,L,K) = &
State_Diag%TomasNUCRATEnumber(I,J,L,K) + (NK(K) - NKD(K)) / &
DTCHEM / BOXMASS ! no./kg air/sec
! print*,'Values for nuc',BOXMASS,DTCHEM,NK(K),NKD(K),K, &
! (NK(K) - NKD(K)), (NK(K) - NKD(K)) / &
! DTCHEM / BOXMASS ! no./kg air/sec
ENDIF
ENDIF
ENDDO
ENDIF

IF ( PTYPE == 11 ) THEN
DO K=1,IBINS
State_Diag%TomasMNFIXezwat1mass(I,J,L,K) = 0.e+0_fp
Expand Down Expand Up @@ -6355,16 +6507,17 @@ SUBROUTINE INIT_TOMAS( Input_Opt, State_Chm, State_Grid, RC )
! (bmy, 1/30/14)
DATA_DIR = TRIM( Input_Opt%CHEM_INPUTS_DIR ) // 'TOMAS_201402/'

! comment out if shut off subgridcoag 14/12/23, bc
! Define subgrid coagulation timescale (win, 10/28/08)
IF ( TRIM(State_Grid%GridRes) == '4.0x5.0' ) THEN
SGCTSCALE = 10.*3600. ! 10 hours
ELSE IF ( TRIM(State_Grid%GridRes) == '2.0x2.5' ) THEN
SGCTSCALE = 5.*3600.
ELSE IF ( TRIM(State_Grid%GridRes) == '0.5x0.625' ) THEN
SGCTSCALE = 1.*3600.
ELSE IF ( TRIM(State_Grid%GridRes) == '0.25x0.3125' ) THEN
SGCTSCALE = 0.5*3600.
ENDIF
!IF ( TRIM(State_Grid%GridRes) == '4.0x5.0' ) THEN
! SGCTSCALE = 10.0_fp*3600.0_fp ! 10 hours
!ELSE IF ( TRIM(State_Grid%GridRes) == '2.0x2.5' ) THEN
! SGCTSCALE = 5.0_fp*3600.0_fp
!ELSE IF ( TRIM(State_Grid%GridRes) == '0.5x0.625' ) THEN
! SGCTSCALE = 1.0_fp*3600.0_fp
!ELSE IF ( TRIM(State_Grid%GridRes) == '0.25x0.3125' ) THEN
! SGCTSCALE = 0.50_fp*3600.0_fp
!ENDIF
msulprizio marked this conversation as resolved.
Show resolved Hide resolved

#if defined(TOMAS40)
Mo = 1.0e-21_fp*2.e+0_fp**(-10)
Expand Down Expand Up @@ -8559,7 +8712,7 @@ SUBROUTINE SUBGRIDCOAG( NDISTINIT, NDIST, MDIST, BOXVOLUME,TEMP, &
if (kcoag(k).gt.0.e+0_fp)then
fracdiaml(k,kk)=kij(k,kk)*ndist(kk)/kcoag(k)
else
fracdiaml(k,kk)=0
fracdiaml(k,kk)=0.e+0_fp
endif
!debug
if(pdbug) print *, kk, fracdiaml(k,kk)
Expand All @@ -8568,7 +8721,9 @@ SUBROUTINE SUBGRIDCOAG( NDISTINIT, NDIST, MDIST, BOXVOLUME,TEMP, &

! determine the number of new particles left after coagulation
do k=1,ibins
ndistfinal(k)=ndistinit(k)*exp(-kcoag(k)*tscale)
!ndistfinal(k)=ndistinit(k)*exp(-kcoag(k)*tscale)
!print*, 'Betty is here tscale ',tscale,k,kcoag(k)
ndistfinal(k)=ndistinit(k)*exp(-1.e+0_fp*kcoag(k)*tscale)
enddo

! determine the mass added to each bin coagulation
Expand Down
2 changes: 1 addition & 1 deletion run/GCClassic/createRunDir.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1045,7 +1045,7 @@ fi

#--------------------------------------------------------------------
# Copy sample restart file to run directory
# Bash functions used here are from ../shared/setupForRestarts.sh
# Bash functions used here are from ./setupForRestarts.sh
#--------------------------------------------------------------------
msulprizio marked this conversation as resolved.
Show resolved Hide resolved

# Parse the download_data.yml file, which returns variable declarations
Expand Down
Loading