From 96e21e875d33e7de4fc2eec5d9173cdce27e020a Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Wed, 3 Mar 2021 11:56:48 -0500 Subject: [PATCH 1/3] Fix for PChem for negative values --- GEOSpchem_GridComp/GEOS_PChemGridComp.F90 | 94 ++++++++++++++++++++--- 1 file changed, 84 insertions(+), 10 deletions(-) diff --git a/GEOSpchem_GridComp/GEOS_PChemGridComp.F90 b/GEOSpchem_GridComp/GEOS_PChemGridComp.F90 index fa2da2fa..adb8ad71 100644 --- a/GEOSpchem_GridComp/GEOS_PChemGridComp.F90 +++ b/GEOSpchem_GridComp/GEOS_PChemGridComp.F90 @@ -1748,6 +1748,18 @@ subroutine UPDATE(NN,NAME,XX) return end if endif + +! if ( ANY(XX < 0.0) ) then +! do l=1,lm +! do j=1,jm +! do i=1,im +! if ( XX(i,j,l) < 0.0 ) then +! print*,'NEG '//trim(NAME)//' AT LEV=', l, XX(i,j,l) +! end if +! enddo +! enddo +! enddo +! endif call MAPL_GetResource(MAPL, TAU,LABEL=trim(NAME)//"_RELAXTIME:", DEFAULT=0.0 ,RC=STATUS) VERIFY_(STATUS) @@ -1777,12 +1789,12 @@ subroutine UPDATE(NN,NAME,XX) do j=1,jm do l=1,nlevs - call MAPL_INTERP( PROD(:,L), LATS(:,J), Prod1(:,L), PCHEM_STATE%LATS) - call MAPL_INTERP( LOSS(:,L), LATS(:,J), Loss1(:,L), PCHEM_STATE%LATS) + call INTERP_NO_EXTRAP( PROD(:,L), LATS(:,J), Prod1(:,L), PCHEM_STATE%LATS) + call INTERP_NO_EXTRAP( LOSS(:,L), LATS(:,J), Loss1(:,L), PCHEM_STATE%LATS) enddo do i=1,im - call MAPL_INTERP( PROD_INT(i,j,:), PL(i,j,:), PROD(i,:), PCHEM_STATE%LEVS) - call MAPL_INTERP( LOSS_INT(i,j,:), PL(i,j,:), LOSS(i,:), PCHEM_STATE%LEVS) + call INTERP_NO_EXTRAP( PROD_INT(i,j,:), PL(i,j,:), PROD(i,:), PCHEM_STATE%LEVS) + call INTERP_NO_EXTRAP( LOSS_INT(i,j,:), PL(i,j,:), LOSS(i,:), PCHEM_STATE%LEVS) enddo end do @@ -1794,10 +1806,10 @@ subroutine UPDATE(NN,NAME,XX) do j=1,jm do l=1,nlevs - call MAPL_INTERP( PROD(:,L), LATS(:,J), Prod1(:,L), PCHEM_STATE%LATS) + call INTERP_NO_EXTRAP( PROD(:,L), LATS(:,J), Prod1(:,L), PCHEM_STATE%LATS) enddo do i=1,im - call MAPL_INTERP( PROD_INT(i,j,:), PL(i,j,:), PROD(i,:), PCHEM_STATE%LEVS) + call INTERP_NO_EXTRAP( PROD_INT(i,j,:), PL(i,j,:), PROD(i,:), PCHEM_STATE%LEVS) enddo end do @@ -1827,6 +1839,29 @@ subroutine UPDATE(NN,NAME,XX) LOSS_INT = (1./TAU) * max( min( (PCRIT -PL)/DELP, 1.0), 0.0) endif +! if ( ANY(PROD_INT < 0.0) ) then +! do l=1,lm +! do j=1,jm +! do i=1,im +! if ( PROD_INT(i,j,l) < 0.0 ) then +! print*,'NEG PROD_INT '//trim(NAME)//' at lev=', l, PROD_INT(i,j,l) +! end if +! enddo +! enddo +! enddo +! endif +! if ( ANY(LOSS_INT < 0.0) ) then +! do l=1,lm +! do j=1,jm +! do i=1,im +! if ( LOSS_INT(i,j,l) < 0.0 ) then +! print*,'NEG LOSS_INT '//trim(NAME)//' at lev=', l, LOSS_INT(i,j,l) +! end if +! enddo +! enddo +! enddo +! endif + PROD_INT = LOSS_INT*PROD_INT XX = (XX + DT*PROD_INT) / (1.0 + DT*LOSS_INT) @@ -1849,6 +1884,18 @@ subroutine UPDATE(NN,NAME,XX) if(associated(H2O_TEND)) H2O_TEND = (PROD_INT - LOSS_INT*XX) end if +! if ( ANY(XX < 0.0) ) then +! do l=1,lm +! do j=1,jm +! do i=1,im +! if ( XX(i,j,l) < 0.0 ) then +! print*,'NEG '//trim(NAME)//' at lev=', l, XX(i,j,l) +! end if +! enddo +! enddo +! enddo +! endif + return end subroutine UPDATE @@ -1885,12 +1932,12 @@ subroutine UPDATE_H2O_PL(NAME,XX) do j=1,jm do l=1,nlevs - call MAPL_INTERP( PROD(:,L), LATS(:,J), Prod1(:,L), PCHEM_STATE%LATS) - call MAPL_INTERP( LOSS(:,L), LATS(:,J), Loss1(:,L), PCHEM_STATE%LATS) + call INTERP_NO_EXTRAP( PROD(:,L), LATS(:,J), Prod1(:,L), PCHEM_STATE%LATS) + call INTERP_NO_EXTRAP( LOSS(:,L), LATS(:,J), Loss1(:,L), PCHEM_STATE%LATS) enddo do i=1,im - call MAPL_INTERP( PROD_INT(i,j,:), PL(i,j,:), PROD(i,:), PCHEM_STATE%LEVS) - call MAPL_INTERP( LOSS_INT(i,j,:), PL(i,j,:), LOSS(i,:), PCHEM_STATE%LEVS) + call INTERP_NO_EXTRAP( PROD_INT(i,j,:), PL(i,j,:), PROD(i,:), PCHEM_STATE%LEVS) + call INTERP_NO_EXTRAP( LOSS_INT(i,j,:), PL(i,j,:), LOSS(i,:), PCHEM_STATE%LEVS) enddo end do @@ -2077,5 +2124,32 @@ subroutine AINC_UPDATE ( GC, IMPORT, EXPORT, CLOCK, RC ) end subroutine AINC_UPDATE + subroutine INTERP_NO_EXTRAP( OY, OX, IY, IX ) + +! A variation of INTERP_LIN_1111_1 from MAPL_InterpMod.F90 +! in which the interpolation is not allowed to extrapolate. +! ASSUMPTION: Values in IX are INCREASING (MAPL assumes this too) + + real, intent(OUT) :: OY(:) + real, intent(IN ) :: OX(:) + real, intent(IN ) :: IY(:) + real, intent(IN ) :: IX(:) + + integer max_index + + max_index = size(IX) + + call MAPL_INTERP( OY, OX, IY, IX ) + + where ( OX < IX(1) ) + OY = IY(1) + end where + + where ( OX > IX(max_index) ) + OY = IY(max_index) + end where + + return + end subroutine INTERP_NO_EXTRAP end module GEOS_PChemGridCompMod From 7f6818569b9d354c6e31d4dfc2f407b8893d8530 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Wed, 3 Mar 2021 12:52:27 -0500 Subject: [PATCH 2/3] Remove commented out debugging code --- GEOSpchem_GridComp/GEOS_PChemGridComp.F90 | 47 ----------------------- 1 file changed, 47 deletions(-) diff --git a/GEOSpchem_GridComp/GEOS_PChemGridComp.F90 b/GEOSpchem_GridComp/GEOS_PChemGridComp.F90 index adb8ad71..9a409ffe 100644 --- a/GEOSpchem_GridComp/GEOS_PChemGridComp.F90 +++ b/GEOSpchem_GridComp/GEOS_PChemGridComp.F90 @@ -1749,18 +1749,6 @@ subroutine UPDATE(NN,NAME,XX) end if endif -! if ( ANY(XX < 0.0) ) then -! do l=1,lm -! do j=1,jm -! do i=1,im -! if ( XX(i,j,l) < 0.0 ) then -! print*,'NEG '//trim(NAME)//' AT LEV=', l, XX(i,j,l) -! end if -! enddo -! enddo -! enddo -! endif - call MAPL_GetResource(MAPL, TAU,LABEL=trim(NAME)//"_RELAXTIME:", DEFAULT=0.0 ,RC=STATUS) VERIFY_(STATUS) @@ -1839,29 +1827,6 @@ subroutine UPDATE(NN,NAME,XX) LOSS_INT = (1./TAU) * max( min( (PCRIT -PL)/DELP, 1.0), 0.0) endif -! if ( ANY(PROD_INT < 0.0) ) then -! do l=1,lm -! do j=1,jm -! do i=1,im -! if ( PROD_INT(i,j,l) < 0.0 ) then -! print*,'NEG PROD_INT '//trim(NAME)//' at lev=', l, PROD_INT(i,j,l) -! end if -! enddo -! enddo -! enddo -! endif -! if ( ANY(LOSS_INT < 0.0) ) then -! do l=1,lm -! do j=1,jm -! do i=1,im -! if ( LOSS_INT(i,j,l) < 0.0 ) then -! print*,'NEG LOSS_INT '//trim(NAME)//' at lev=', l, LOSS_INT(i,j,l) -! end if -! enddo -! enddo -! enddo -! endif - PROD_INT = LOSS_INT*PROD_INT XX = (XX + DT*PROD_INT) / (1.0 + DT*LOSS_INT) @@ -1884,18 +1849,6 @@ subroutine UPDATE(NN,NAME,XX) if(associated(H2O_TEND)) H2O_TEND = (PROD_INT - LOSS_INT*XX) end if -! if ( ANY(XX < 0.0) ) then -! do l=1,lm -! do j=1,jm -! do i=1,im -! if ( XX(i,j,l) < 0.0 ) then -! print*,'NEG '//trim(NAME)//' at lev=', l, XX(i,j,l) -! end if -! enddo -! enddo -! enddo -! endif - return end subroutine UPDATE From 9b2a16f97a857d97a34bb1f8567b1ede808ab945 Mon Sep 17 00:00:00 2001 From: Matthew Thompson Date: Tue, 6 Apr 2021 14:15:46 -0400 Subject: [PATCH 3/3] Update CI image to 6.1.0 --- .circleci/config.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.circleci/config.yml b/.circleci/config.yml index f9781817..c2dc41bd 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -3,7 +3,7 @@ version: 2.1 executors: gcc-build-env: docker: - - image: gmao/ubuntu20-geos-env-mkl:v6.0.27-openmpi_4.0.5-gcc_10.2.0 + - image: gmao/ubuntu20-geos-env-mkl:v6.1.0-openmpi_4.0.5-gcc_10.2.0 auth: username: $DOCKERHUB_USER password: $DOCKERHUB_AUTH_TOKEN