Back to home page

MITgcm

 
 

    


File indexing completed on 2018-03-02 18:40:43 UTC

view on githubraw file Latest commit add29e06 on 2018-01-31 20:35:05 UTC
a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
e337e4ca8c Andr*0002        subroutine update_chemistry_exports (myTime, myIter, myThid)
                0003 c----------------------------------------------------------------------
                0004 c  Subroutine update_chemistry_exports - 'Wrapper' routine to update
3daafce20b Jean*0005 c        the fields related to the earth chemistry that are needed
e337e4ca8c Andr*0006 c        by fizhi.
                0007 c        Also: Set up "bi, bj loop" and some timers and clocks here.
                0008 c
                0009 c Call:  interp_chemistry
                0010 c-----------------------------------------------------------------------
                0011        implicit none
                0012 #include "SIZE.h"
                0013 #include "fizhi_SIZE.h"
613fa3996d Andr*0014 #include "fizhi_land_SIZE.h"
e337e4ca8c Andr*0015 #include "GRID.h"
                0016 #include "DYNVARS.h"
f4a0368053 Andr*0017 #include "fizhi_chemistry_coms.h"
86214b2dce Andr*0018 #include "fizhi_coms.h"
e337e4ca8c Andr*0019 #include "gridalt_mapping.h"
                0020 #include "EEPARAMS.h"
86214b2dce Andr*0021 #include "chronos.h"
e337e4ca8c Andr*0022 
3768927683 Andr*0023       integer myIter, myThid
                0024       _RL myTime
e337e4ca8c Andr*0025 
                0026 c pe on physics grid refers to bottom edge
a456aa407c Andr*0027       _RL pephy(1-Olx:sNx+Olx,1-Oly:sNy+Oly,Nrphys+1,nSx,nSy)
                0028       _RL pphy(sNx,sNy,Nrphys,nSx,nSy)
                0029       _RL oz1(nlatsoz,nlevsoz), strq1(nlatsq,nlevsq)
                0030       _RL waterin(sNx,sNy,Nrphys), xlat(sNx,sNy)
86214b2dce Andr*0031 
d0b11e35fb Andr*0032       integer i, j, L, LL, bi, bj
9524fe64b5 Andr*0033       integer im1, im2, jm1, jm2
86214b2dce Andr*0034       integer nhms1,nymd1,nhms2,nymd2,imns,ipls
a456aa407c Andr*0035       _RL facm, facp
2a3ae9099b Andr*0036       logical alarm
                0037       external alarm
86214b2dce Andr*0038 
9524fe64b5 Andr*0039       im1 = 1
                0040       im2 = sNx
                0041       jm1 = 1
                0042       jm2 = sNy
86214b2dce Andr*0043 
                0044       if( alarm('radsw').or.alarm('radlw') ) then
                0045 
e337e4ca8c Andr*0046        do bj = myByLo(myThid), myByHi(myThid)
                0047        do bi = myBxLo(myThid), myBxHi(myThid)
                0048 
d0b11e35fb Andr*0049 c  Construct the physics grid pressures - count pephy levels top down
                0050 c                                         (even though dpphy counted bottom up)
e337e4ca8c Andr*0051         do j = 1,sNy
                0052         do i = 1,sNx
00f44e1146 Andr*0053          pephy(i,j,Nrphys+1,bi,bj)=(Ro_surf(i,j,bi,bj)+etaH(i,j,bi,bj))
e337e4ca8c Andr*0054          do L = 2,Nrphys+1
d0b11e35fb Andr*0055          LL = Nrphys+2-L
                0056          pephy(i,j,LL,bi,bj)=pephy(i,j,LL+1,bi,bj)-dpphys(i,j,L-1,bi,bj)
e337e4ca8c Andr*0057          enddo
                0058         enddo
                0059         enddo
86214b2dce Andr*0060         do j = 1,sNy
                0061         do i = 1,sNx
                0062          do L = 1,Nrphys
9524fe64b5 Andr*0063           pphy(i,j,L,bi,bj)=(pephy(i,j,L+1,bi,bj)+pephy(i,j,L,bi,bj))
                0064      .                                              /200.
86214b2dce Andr*0065          enddo
                0066         enddo
                0067         enddo
                0068 
                0069         do j = 1,sNy
                0070         do i = 1,sNx
                0071          xlat(i,j) = yC(i,j,bi,bj)
                0072          do L = 1,Nrphys
                0073           waterin(i,j,L) = sphy(i,j,L,bi,bj)
                0074          enddo
                0075         enddo
                0076         enddo
                0077 
                0078         call time_bound(nymd,nhms,nymd1,nhms1,nymd2,nhms2,imns,ipls)
                0079         call interp_time(nymd,nhms,nymd1,nhms1,nymd2,nhms2,facm,facp)
                0080 
                0081         do L = 1,nlevsoz
                0082         do j = 1,nlatsoz
                0083          oz1(j,L) = ozone(j,L,imns)*facm + ozone(j,L,ipls)*facp
                0084         enddo
                0085         enddo
0889f02121 Jean*0086 
86214b2dce Andr*0087         do L = 1,nlevsq
                0088         do j = 1,nlatsq
                0089          strq1(j,L) = stratq(j,L,imns)*facm + stratq(j,L,ipls)*facp
                0090         enddo
                0091         enddo
e337e4ca8c Andr*0092 
86214b2dce Andr*0093         call interp_chemistry(strq1,nlevsq,nlatsq,levsq,latsq,
33b0afdcb4 Jean*0094      .   oz1,nlevsoz,nlatsoz,levsoz,latsoz,
                0095      .   waterin,pphy(1,1,1,bi,bj),xlat,
86214b2dce Andr*0096      .   im2,jm2,Nrphys,nSx,nSy,bi,bj,o3,qstr)
e337e4ca8c Andr*0097 
                0098        enddo
                0099        enddo
                0100 
86214b2dce Andr*0101       endif
                0102 
                0103       return
                0104       end
63416ca6a5 Andr*0105 
                0106       subroutine interp_chemistry (stratq,nwatlevs,nwatlats,watlevs,
613fa3996d Andr*0107      . watlats,ozone,nozlevs,nozlats,ozlevs,ozlats,
                0108      . qz,plz,xlat,im,jm,lm,nSx,nSy,bi,bj,ozrad,qzrad)
63416ca6a5 Andr*0109 
                0110       implicit none
                0111 
                0112 c Input Variables
                0113 c ---------------
613fa3996d Andr*0114       integer nwatlevs,nwatlats,nozlevs,nozlats,nSx,nSy,bi,bj
a456aa407c Andr*0115       _RL stratq(nwatlats,nwatlevs),ozone(nozlats,nozlevs)
613fa3996d Andr*0116       _RL watlevs(nwatlevs),watlats(nwatlats)
                0117       _RL ozlevs(nozlevs),ozlats(nozlats)
63416ca6a5 Andr*0118       integer im,jm,lm
a456aa407c Andr*0119       _RL qz(im,jm,lm),plz(im,jm,lm)
                0120       _RL xlat(im,jm)
613fa3996d Andr*0121       _RL ozrad(im,jm,lm,nSx,nSy)
                0122       _RL qzrad(im,jm,lm,nSx,nSy)
63416ca6a5 Andr*0123 
                0124 C **********************************************************************
                0125 C ****           Get Ozone and Stratospheric Moisture Data          ****
                0126 C **********************************************************************
                0127 
                0128       call interp_qz (stratq,nwatlevs,nwatlats,watlevs,watlats,im*jm,
9524fe64b5 Andr*0129      .                         bi,bj, xlat,lm,plz,qz,qzrad(1,1,1,bi,bj))
613fa3996d Andr*0130       call interp_oz (ozone ,nozlevs,nozlats,ozlevs,ozlats,im*jm,
9524fe64b5 Andr*0131      .                         bi,bj, xlat,lm,plz,ozrad(1,1,1,bi,bj))
                0132 
63416ca6a5 Andr*0133       return
                0134       end
                0135 
                0136       subroutine interp_qz(stratq,nwatlevs,nwatlats,watlevs,watlats,
9524fe64b5 Andr*0137      .                         irun,bi,bj,xlat,nlevs,pres,qz_in,qz_out )
63416ca6a5 Andr*0138 C***********************************************************************
                0139 C  Purpose
                0140 C     To Interpolate Chemistry Moisture from Chemistry Grid to Physics Grid
                0141 C
                0142 C  INPUT Argument Description
                0143 C     stratq .... Climatological SAGE Stratospheric Moisture
                0144 C     irun ...... Number of Columns to be filled
                0145 C     xlat ...... Latitude in Degrees
                0146 C     nlevs ..... Vertical   Dimension
                0147 C     pres ...... PRES (IM,JM,nlevs) Three-dimensional array of pressures
                0148 C     qz_in ..... Model Moisture (kg/kg mass mixing radtio)
0889f02121 Jean*0149 C     qz_out .... Combination of Chemistry Moisture and Model Moisture
613fa3996d Andr*0150 C                 (kg/kg mass mixing ratio)
63416ca6a5 Andr*0151 C
                0152 C***********************************************************************
                0153 
0889f02121 Jean*0154       implicit none
63416ca6a5 Andr*0155       integer nwatlevs,nwatlats
9524fe64b5 Andr*0156       integer bi,bj
a456aa407c Andr*0157       _RL stratq ( nwatlats,nwatlevs )
613fa3996d Andr*0158       _RL watlats (nwatlats)
                0159       _RL watlevs (nwatlevs)
63416ca6a5 Andr*0160 
                0161       integer irun,nlevs
a456aa407c Andr*0162       _RL xlat  (irun)
                0163       _RL pres  (irun,nlevs)
                0164       _RL qz_in (irun,nlevs)
613fa3996d Andr*0165       _RL qz_out(irun,nlevs)
63416ca6a5 Andr*0166 
                0167 c Local Variables
                0168 c ---------------
                0169       integer     pqu,pql,dpq
                0170       parameter ( pqu = 100.    )
                0171       parameter ( pql = 300.    )
                0172       parameter ( dpq = pql-pqu )
                0173 
                0174       integer i,k,L1,L2,LM,LP
a456aa407c Andr*0175       _RL h2o_time_lat (irun,nwatlevs)
                0176       _RL       qz_clim(irun,nlevs)
63416ca6a5 Andr*0177 
a456aa407c Andr*0178       _RL  qpr1(irun), qpr2(irun), slope(irun)
                0179       _RL   pr1(irun),  pr2(irun)
63416ca6a5 Andr*0180 
                0181       integer  jlat,jlatm,jlatp
                0182 
                0183 C **********************************************************************
                0184 C ****         Interpolate Moisture data to model latitudes          ***
                0185 C **********************************************************************
0889f02121 Jean*0186 
63416ca6a5 Andr*0187       DO 32 k = 1, nwatlevs
                0188         DO 34   i = 1,irun
                0189 
                0190         DO 36 jlat = 1, nwatlats
                0191            IF( watlats(jlat).gt.xlat(i) ) THEN
                0192               IF( jlat.EQ.1 ) THEN
                0193                   jlatm    = 1
                0194                   jlatp    = 1
                0195                   slope(i) = 0
                0196                     ELSE
                0197                   jlatm    = jlat -1
                0198                   jlatp    = jlat
                0199                   slope(i) = ( xlat(i)        -watlats(jlat-1) )
                0200      .                     / ( watlats(jlat)-watlats(jlat-1) )
                0201               ENDIF
                0202               GOTO 37
                0203            ENDIF
                0204    36   CONTINUE
                0205         jlatm    = nwatlats
                0206         jlatp    = nwatlats
                0207         slope(i) =  1
                0208    37   CONTINUE
                0209         QPR1(i) = stratq(jlatm,k)
                0210         QPR2(i) = stratq(jlatp,k)
                0211    34   CONTINUE
                0212 
                0213         do  i = 1,irun
                0214         h2o_time_lat(i,k) = qpr1(i) + slope(i)*(qpr2(i)-qpr1(i))
                0215         enddo
                0216 
                0217    32 CONTINUE
                0218 
                0219 C **********************************************************************
                0220 C ****     Interpolate Latitude Moisture data to model pressures     ***
                0221 C **********************************************************************
0889f02121 Jean*0222 
63416ca6a5 Andr*0223       DO 40 L2 = 1,nlevs
                0224 
                0225         DO 44 i= 1, irun
                0226         DO 46 L1 = 1,nwatlevs
                0227            IF( watlevs(L1).GT.pres(i,L2) ) THEN
                0228              IF( L1.EQ.1 ) THEN
                0229                  LM = 1
                0230                  LP = 2
                0231                ELSE
                0232                  LM = L1-1
                0233                  LP = L1
                0234              ENDIF
                0235              GOTO 47
                0236            ENDIF
                0237    46   CONTINUE
                0238         LM = nwatlevs-1
                0239         LP = nwatlevs
                0240    47   CONTINUE
                0241          PR1(i) =     watlevs (LM)
                0242          PR2(i) =     watlevs (LP)
                0243         QPR1(i) = h2o_time_lat(i,LM)
                0244         QPR2(i) = h2o_time_lat(i,LP)
                0245    44   CONTINUE
                0246 
                0247       do i= 1, irun
                0248            slope(i) =(QPR1(i)-QPR2(i)) / (PR1(i)-PR2(i))
                0249       qz_clim(i,L2) = QPR2(i) + (pres(i,L2)-PR2(i))*SLOPE(i)
                0250       enddo
                0251 
                0252    40 CONTINUE
                0253 
                0254 c
                0255 c ... Above 100 mb, using climatological  water data set ...................
                0256 c ... Below 300 mb, using model predicted water data set ...................
                0257 c ... In between, using linear interpolation ...............................
                0258 c
                0259       do k= 1, nlevs
                0260       do i= 1, irun
                0261            if( pres(i,k).ge.pqu  .and. pres(i,k).le. pql) then
                0262              qz_out(i,k) = qz_clim(i,k)+(qz_in(i,k)-
                0263      1                     qz_clim(i,k))*(pres(i,k)-pqu)/dpq
                0264       else if( pres(i,k) .gt. pql ) then
                0265              qz_out(i,k) = qz_in  (i,k)
                0266       else
                0267              qz_out(i,k) = qz_clim(i,k)
                0268            endif
                0269       enddo
                0270       enddo
                0271 
                0272       return
                0273       end
                0274 
613fa3996d Andr*0275       subroutine interp_oz (ozone,nozlevs,nozlats,ozlevs,ozlats,
9524fe64b5 Andr*0276      .                            irun,bi,bj,xlat,nlevs,plevs,ozrad)
63416ca6a5 Andr*0277 C***********************************************************************
                0278 C  Purpose
                0279 C     To Interpolate Chemistry Ozone from Chemistry Grid to Physics Grid
                0280 C
                0281 C  INPUT Argument Description
                0282 C     ozone ..... Climatological Ozone
                0283 C     chemistry .. Chemistry State Data Structure
                0284 C     irun ....... Number of Columns to be filled
                0285 C     xlat ....... Latitude in Degrees
                0286 C     nlevs ...... Vertical   Dimension
                0287 C     pres ....... Three-dimensional array of pressures
                0288 C     ozrad ...... Ozone on Physics Grid (kg/kg mass mixing radtio)
                0289 C
                0290 C***********************************************************************
                0291       implicit none
613fa3996d Andr*0292       integer nozlevs,nozlats,irun,nlevs
9524fe64b5 Andr*0293       integer bi,bj
a456aa407c Andr*0294       _RL ozone(nozlats,nozlevs)
                0295       _RL xlat(irun)
                0296       _RL plevs(irun,nlevs)
613fa3996d Andr*0297       _RL ozrad(irun,nlevs)
                0298       _RL ozlevs(nozlevs),ozlats(nozlats)
63416ca6a5 Andr*0299 
                0300 c Local Variables
                0301 c ---------------
a456aa407c Andr*0302       _RL zero,one,o3min,voltomas
63416ca6a5 Andr*0303       PARAMETER ( ZERO     = 0.0 )
                0304       PARAMETER ( ONE      = 1.0 )
                0305       PARAMETER ( O3MIN    = 1.0E-10  )
                0306       PARAMETER ( VOLTOMAS = 1.655E-6 )
                0307 
                0308       integer  i,k,L1,L2,LM,LP
                0309       integer  jlat,jlatm,jlatp
a456aa407c Andr*0310       _RL  O3INT1(IRUN,nozlevs)
                0311       _RL    QPR1(IRUN), QPR2(IRUN), SLOPE(IRUN)
                0312       _RL     PR1(IRUN),  PR2(IRUN)
63416ca6a5 Andr*0313 
                0314 C **********************************************************************
                0315 C ****           INTERPOLATE ozone data to model latitudes           ***
                0316 C **********************************************************************
                0317 
613fa3996d Andr*0318       DO 32 K=1,nozlevs
63416ca6a5 Andr*0319       DO 34 I=1,IRUN
                0320 
613fa3996d Andr*0321       DO 36 jlat = 1,nozlats
                0322       IF( ozlats(jlat).gt.xlat(i) ) THEN
63416ca6a5 Andr*0323       IF( jlat.EQ.1 ) THEN
                0324       jlatm    = 1
                0325       jlatp    = 1
                0326       slope(i) = zero
                0327         ELSE
                0328       jlatm    = jlat-1
                0329       jlatp    = jlat
613fa3996d Andr*0330       slope(i) = ( XLAT(I)        -ozlats(jlat-1) )
                0331      .         / ( ozlats(jlat)-ozlats(jlat-1) )
63416ca6a5 Andr*0332       ENDIF
                0333       GOTO 37
                0334       ENDIF
                0335    36 CONTINUE
613fa3996d Andr*0336       jlatm    = nozlats
                0337       jlatp    = nozlats
63416ca6a5 Andr*0338       slope(i) = one
                0339    37 CONTINUE
                0340       QPR1(I) = ozone(jlatm,k)
                0341       QPR2(I) = ozone(jlatp,k)
                0342    34 CONTINUE
                0343 
                0344       DO 38 I=1,IRUN
                0345       o3int1(i,k) = qpr1(i) + slope(i)*( qpr2(i)-qpr1(i) )
                0346    38 CONTINUE
                0347 
                0348    32 CONTINUE
                0349 
                0350 C **********************************************************************
                0351 C ****     INTERPOLATE latitude ozone data to model pressures        ***
                0352 C **********************************************************************
                0353 
                0354       DO 40 L2 = 1,NLEVS
                0355 
                0356       DO 44 I  = 1,IRUN
613fa3996d Andr*0357       DO 46 L1 = 1,nozlevs
                0358       IF( ozlevs(L1).GT.PLEVS(I,L2) ) THEN
63416ca6a5 Andr*0359       IF( L1.EQ.1 ) THEN
                0360           LM = 1
                0361           LP = 2
                0362         ELSE
                0363           LM = L1-1
                0364           LP = L1
                0365       ENDIF
                0366       GOTO 47
                0367       ENDIF
                0368    46 CONTINUE
613fa3996d Andr*0369             LM = nozlevs-1
                0370             LP = nozlevs
63416ca6a5 Andr*0371    47 CONTINUE
613fa3996d Andr*0372        PR1(I) = ozlevs (LM)
                0373        PR2(I) = ozlevs (LP)
63416ca6a5 Andr*0374       QPR1(I) =   O3INT1(I,LM)
                0375       QPR2(I) =   O3INT1(I,LP)
                0376    44 CONTINUE
                0377 
                0378       DO 48 I=1,IRUN
                0379          SLOPE(I) = ( QPR1(I)-QPR2(I) )
                0380      .            / (  PR1(I)- PR2(I) )
                0381       ozrad(I,L2) =   QPR2(I) + ( PLEVS(I,L2)-PR2(I) )*SLOPE(I)
                0382 
                0383       if( ozrad(i,l2).lt.o3min ) then
                0384           ozrad(i,l2) =  o3min
                0385       endif
                0386 
                0387    48 CONTINUE
                0388    40 CONTINUE
                0389 
                0390 C **********************************************************************
                0391 C ****     CONVERT FROM VOLUME MIXING RATIO TO MASS MIXING RATIO     ***
                0392 C **********************************************************************
                0393 
0889f02121 Jean*0394       DO 60 L2=1,NLEVS
                0395       DO 60 I=1,IRUN
                0396 c     DO 60 I=1,IRUN*NLEVS
                0397 c     ozrad (I,1) = ozrad(I,1) * VOLTOMAS
                0398       ozrad (I,L2) = ozrad(I,L2) * VOLTOMAS
63416ca6a5 Andr*0399   60  CONTINUE
                0400 
                0401       RETURN
                0402       END
                0403