Merge pull request #125 from GEOS-ESM/develop

Fix PCHEM interpolation (dev->main)
GEOS-ESM · Apr 15, 2021 · 869541f · 869541f
2 parents 7b78f24 + f7a3d50
commit 869541f
Show file tree

Hide file tree

Showing 2 changed files with 39 additions and 12 deletions.
diff --git 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

diff --git a/GEOSpchem_GridComp/GEOS_PChemGridComp.F90 b/GEOSpchem_GridComp/GEOS_PChemGridComp.F90
@@ -1748,7 +1748,7 @@ subroutine UPDATE(NN,NAME,XX)
           return
        end if
     endif
- 
+
     call MAPL_GetResource(MAPL,   TAU,LABEL=trim(NAME)//"_RELAXTIME:", DEFAULT=0.0 ,RC=STATUS)
     VERIFY_(STATUS)
 
@@ -1777,12 +1777,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 +1794,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
 
@@ -1885,12 +1885,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 +2077,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