C $Header: /u/gcmpack/MITgcm/pkg/fizhi/update_chemistry_exports.F,v 1.14 2010/03/16 00:19:33 jmc Exp $
C $Name: $
#include "FIZHI_OPTIONS.h"
subroutine UPDATE_CHEMISTRY_EXPORTS (myTime, myIter, myThid)
c----------------------------------------------------------------------
c Subroutine update_chemistry_exports - 'Wrapper' routine to update
c the fields related to the earth chemistry that are needed
c by fizhi.
c Also: Set up "bi, bj loop" and some timers and clocks here.
c
c Call: interp_chemistry
c-----------------------------------------------------------------------
implicit none
#include "SIZE.h"
#include "fizhi_SIZE.h"
#include "fizhi_land_SIZE.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "fizhi_chemistry_coms.h"
#include "fizhi_coms.h"
#include "gridalt_mapping.h"
#include "EEPARAMS.h"
#include "chronos.h"
integer myIter, myThid
_RL myTime
c pe on physics grid refers to bottom edge
_RL pephy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys+1,nSx,nSy)
_RL pphy(sNx,sNy,Nrphys,nSx,nSy)
_RL oz1(nlatsoz,nlevsoz), strq1(nlatsq,nlevsq)
_RL waterin(sNx,sNy,Nrphys), xlat(sNx,sNy)
integer i, j, L, LL, bi, bj
integer im1, im2, jm1, jm2
integer nhms1,nymd1,nhms2,nymd2,imns,ipls
_RL facm, facp
logical alarm
external
im1 = 1
im2 = sNx
jm1 = 1
jm2 = sNy
if( alarm('radsw').or.alarm('radlw') ) then
do bj = myByLo(myThid), myByHi(myThid)
do bi = myBxLo(myThid), myBxHi(myThid)
c Construct the physics grid pressures - count pephy levels top down
c (even though dpphy counted bottom up)
do j = 1,sNy
do i = 1,sNx
pephy(i,j,Nrphys+1,bi,bj)=(Ro_surf(i,j,bi,bj)+etaH(i,j,bi,bj))
do L = 2,Nrphys+1
LL = Nrphys+2-L
pephy(i,j,LL,bi,bj)=pephy(i,j,LL+1,bi,bj)-dpphys(i,j,L-1,bi,bj)
enddo
enddo
enddo
do j = 1,sNy
do i = 1,sNx
do L = 1,Nrphys
pphy(i,j,L,bi,bj)=(pephy(i,j,L+1,bi,bj)+pephy(i,j,L,bi,bj))
. /200.
enddo
enddo
enddo
do j = 1,sNy
do i = 1,sNx
xlat(i,j) = yC(i,j,bi,bj)
do L = 1,Nrphys
waterin(i,j,L) = sphy(i,j,L,bi,bj)
enddo
enddo
enddo
call TIME_BOUND(nymd,nhms,nymd1,nhms1,nymd2,nhms2,imns,ipls)
call INTERP_TIME(nymd,nhms,nymd1,nhms1,nymd2,nhms2,facm,facp)
do L = 1,nlevsoz
do j = 1,nlatsoz
oz1(j,L) = ozone(j,L,imns)*facm + ozone(j,L,ipls)*facp
enddo
enddo
do L = 1,nlevsq
do j = 1,nlatsq
strq1(j,L) = stratq(j,L,imns)*facm + stratq(j,L,ipls)*facp
enddo
enddo
call INTERP_CHEMISTRY(strq1,nlevsq,nlatsq,levsq,latsq,
. oz1,nlevsoz,nlatsoz,levsoz,latsoz,
. waterin,pphy(1,1,1,bi,bj),xlat,
. im2,jm2,Nrphys,nSx,nSy,bi,bj,o3,qstr)
enddo
enddo
endif
return
end
subroutine INTERP_CHEMISTRY (stratq,nwatlevs,nwatlats,watlevs,
. watlats,ozone,nozlevs,nozlats,ozlevs,ozlats,
. qz,plz,xlat,im,jm,lm,nSx,nSy,bi,bj,ozrad,qzrad)
implicit none
c Input Variables
c ---------------
integer nwatlevs,nwatlats,nozlevs,nozlats,nSx,nSy,bi,bj
_RL stratq(nwatlats,nwatlevs),ozone(nozlats,nozlevs)
_RL watlevs(nwatlevs),watlats(nwatlats)
_RL ozlevs(nozlevs),ozlats(nozlats)
integer im,jm,lm
_RL qz(im,jm,lm),plz(im,jm,lm)
_RL xlat(im,jm)
_RL ozrad(im,jm,lm,nSx,nSy)
_RL qzrad(im,jm,lm,nSx,nSy)
C **********************************************************************
C **** Get Ozone and Stratospheric Moisture Data ****
C **********************************************************************
call INTERP_QZ (stratq,nwatlevs,nwatlats,watlevs,watlats,im*jm,
. bi,bj, xlat,lm,plz,qz,qzrad(1,1,1,bi,bj))
call INTERP_OZ (ozone ,nozlevs,nozlats,ozlevs,ozlats,im*jm,
. bi,bj, xlat,lm,plz,ozrad(1,1,1,bi,bj))
return
end
subroutine INTERP_QZ(stratq,nwatlevs,nwatlats,watlevs,watlats,
. irun,bi,bj,xlat,nlevs,pres,qz_in,qz_out )
C***********************************************************************
C Purpose
C To Interpolate Chemistry Moisture from Chemistry Grid to Physics Grid
C
C INPUT Argument Description
C stratq .... Climatological SAGE Stratospheric Moisture
C irun ...... Number of Columns to be filled
C xlat ...... Latitude in Degrees
C nlevs ..... Vertical Dimension
C pres ...... PRES (IM,JM,nlevs) Three-dimensional array of pressures
C qz_in ..... Model Moisture (kg/kg mass mixing radtio)
C qz_out .... Combination of Chemistry Moisture and Model Moisture
C (kg/kg mass mixing ratio)
C
C***********************************************************************
implicit none
integer nwatlevs,nwatlats
integer bi,bj
_RL stratq ( nwatlats,nwatlevs )
_RL watlats (nwatlats)
_RL watlevs (nwatlevs)
integer irun,nlevs
_RL xlat (irun)
_RL pres (irun,nlevs)
_RL qz_in (irun,nlevs)
_RL qz_out(irun,nlevs)
c Local Variables
c ---------------
integer pqu,pql,dpq
parameter ( pqu = 100. )
parameter ( pql = 300. )
parameter ( dpq = pql-pqu )
integer i,k,L1,L2,LM,LP
_RL h2o_time_lat (irun,nwatlevs)
_RL qz_clim(irun,nlevs)
_RL qpr1(irun), qpr2(irun), slope(irun)
_RL pr1(irun), pr2(irun)
integer jlat,jlatm,jlatp
C **********************************************************************
C **** Interpolate Moisture data to model latitudes ***
C **********************************************************************
DO 32 k = 1, nwatlevs
DO 34 i = 1,irun
DO 36 jlat = 1, nwatlats
IF( watlats(jlat).gt.xlat(i) ) THEN
IF( jlat.EQ.1 ) THEN
jlatm = 1
jlatp = 1
slope(i) = 0
ELSE
jlatm = jlat -1
jlatp = jlat
slope(i) = ( xlat(i) -watlats(jlat-1) )
. / ( watlats(jlat)-watlats(jlat-1) )
ENDIF
GOTO 37
ENDIF
36 CONTINUE
jlatm = nwatlats
jlatp = nwatlats
slope(i) = 1
37 CONTINUE
QPR1(i) = stratq(jlatm,k)
QPR2(i) = stratq(jlatp,k)
34 CONTINUE
do i = 1,irun
h2o_time_lat(i,k) = qpr1(i) + slope(i)*(qpr2(i)-qpr1(i))
enddo
32 CONTINUE
C **********************************************************************
C **** Interpolate Latitude Moisture data to model pressures ***
C **********************************************************************
DO 40 L2 = 1,nlevs
DO 44 i= 1, irun
DO 46 L1 = 1,nwatlevs
IF( watlevs(L1).GT.pres(i,L2) ) THEN
IF( L1.EQ.1 ) THEN
LM = 1
LP = 2
ELSE
LM = L1-1
LP = L1
ENDIF
GOTO 47
ENDIF
46 CONTINUE
LM = nwatlevs-1
LP = nwatlevs
47 CONTINUE
PR1(i) = watlevs (LM)
PR2(i) = watlevs (LP)
QPR1(i) = h2o_time_lat(i,LM)
QPR2(i) = h2o_time_lat(i,LP)
44 CONTINUE
do i= 1, irun
slope(i) =(QPR1(i)-QPR2(i)) / (PR1(i)-PR2(i))
qz_clim(i,L2) = QPR2(i) + (pres(i,L2)-PR2(i))*SLOPE(i)
enddo
40 CONTINUE
c
c ... Above 100 mb, using climatological water data set ...................
c ... Below 300 mb, using model predicted water data set ...................
c ... In between, using linear interpolation ...............................
c
do k= 1, nlevs
do i= 1, irun
if( pres(i,k).ge.pqu .and. pres(i,k).le. pql) then
qz_out(i,k) = qz_clim(i,k)+(qz_in(i,k)-
1 qz_clim(i,k))*(pres(i,k)-pqu)/dpq
else if( pres(i,k) .gt. pql ) then
qz_out(i,k) = qz_in (i,k)
else
qz_out(i,k) = qz_clim(i,k)
endif
enddo
enddo
return
end
subroutine INTERP_OZ (ozone,nozlevs,nozlats,ozlevs,ozlats,
. irun,bi,bj,xlat,nlevs,plevs,ozrad)
C***********************************************************************
C Purpose
C To Interpolate Chemistry Ozone from Chemistry Grid to Physics Grid
C
C INPUT Argument Description
C ozone ..... Climatological Ozone
C chemistry .. Chemistry State Data Structure
C irun ....... Number of Columns to be filled
C xlat ....... Latitude in Degrees
C nlevs ...... Vertical Dimension
C pres ....... Three-dimensional array of pressures
C ozrad ...... Ozone on Physics Grid (kg/kg mass mixing radtio)
C
C***********************************************************************
implicit none
integer nozlevs,nozlats,irun,nlevs
integer bi,bj
_RL ozone(nozlats,nozlevs)
_RL xlat(irun)
_RL plevs(irun,nlevs)
_RL ozrad(irun,nlevs)
_RL ozlevs(nozlevs),ozlats(nozlats)
c Local Variables
c ---------------
_RL zero,one,o3min,voltomas
PARAMETER ( ZERO = 0.0 )
PARAMETER ( ONE = 1.0 )
PARAMETER ( O3MIN = 1.0E-10 )
PARAMETER ( VOLTOMAS = 1.655E-6 )
integer i,k,L1,L2,LM,LP
integer jlat,jlatm,jlatp
_RL O3INT1(IRUN,nozlevs)
_RL QPR1(IRUN), QPR2(IRUN), SLOPE(IRUN)
_RL PR1(IRUN), PR2(IRUN)
C **********************************************************************
C **** INTERPOLATE ozone data to model latitudes ***
C **********************************************************************
DO 32 K=1,nozlevs
DO 34 I=1,IRUN
DO 36 jlat = 1,nozlats
IF( ozlats(jlat).gt.xlat(i) ) THEN
IF( jlat.EQ.1 ) THEN
jlatm = 1
jlatp = 1
slope(i) = zero
ELSE
jlatm = jlat-1
jlatp = jlat
slope(i) = ( XLAT(I) -ozlats(jlat-1) )
. / ( ozlats(jlat)-ozlats(jlat-1) )
ENDIF
GOTO 37
ENDIF
36 CONTINUE
jlatm = nozlats
jlatp = nozlats
slope(i) = one
37 CONTINUE
QPR1(I) = ozone(jlatm,k)
QPR2(I) = ozone(jlatp,k)
34 CONTINUE
DO 38 I=1,IRUN
o3int1(i,k) = qpr1(i) + slope(i)*( qpr2(i)-qpr1(i) )
38 CONTINUE
32 CONTINUE
C **********************************************************************
C **** INTERPOLATE latitude ozone data to model pressures ***
C **********************************************************************
DO 40 L2 = 1,NLEVS
DO 44 I = 1,IRUN
DO 46 L1 = 1,nozlevs
IF( ozlevs(L1).GT.PLEVS(I,L2) ) THEN
IF( L1.EQ.1 ) THEN
LM = 1
LP = 2
ELSE
LM = L1-1
LP = L1
ENDIF
GOTO 47
ENDIF
46 CONTINUE
LM = nozlevs-1
LP = nozlevs
47 CONTINUE
PR1(I) = ozlevs (LM)
PR2(I) = ozlevs (LP)
QPR1(I) = O3INT1(I,LM)
QPR2(I) = O3INT1(I,LP)
44 CONTINUE
DO 48 I=1,IRUN
SLOPE(I) = ( QPR1(I)-QPR2(I) )
. / ( PR1(I)- PR2(I) )
ozrad(I,L2) = QPR2(I) + ( PLEVS(I,L2)-PR2(I) )*SLOPE(I)
if( ozrad(i,l2).lt.o3min ) then
ozrad(i,l2) = o3min
endif
48 CONTINUE
40 CONTINUE
C **********************************************************************
C **** CONVERT FROM VOLUME MIXING RATIO TO MASS MIXING RATIO ***
C **********************************************************************
DO 60 L2=1,NLEVS
DO 60 I=1,IRUN
c DO 60 I=1,IRUN*NLEVS
c ozrad (I,1) = ozrad(I,1) * VOLTOMAS
ozrad (I,L2) = ozrad(I,L2) * VOLTOMAS
60 CONTINUE
RETURN
END