pkg/fizhi/fizhi_init_chem.F

a456aa407c Andr*0001 #include "FIZHI_OPTIONS.h"
c915144f5d Andr*0002       subroutine fizhi_init_chem(mythid,nozlats,nozlevs,ntimesoz,
                0003      . ozlats,ozlevs,o3,
                0004      . nwatlats,nwatlevs,ntimesq,watlats,watlevs,water,
9a6b9d7b6d Andr*0005      . Nrphys,pressure,n2o,methane,co2,cfc11,cfc12,cfc22)
328acc481c Andr*0006 C***********************************************************************
                0007 C Subroutine fizhi_init_chem - routine to read in the ozone and upper
                0008 C      atmosphere water vapor, and set the methane, n2o and cfc values
                0009 C
                0010 C INPUT:
                0011 C
                0012 C mythid   - thread number (processor number)
                0013 C
                0014 C OUTPUT:
                0015 C
                0016 C***********************************************************************
                0017       implicit none
                0018       integer mythid,nozlevs,nozlats,nwatlevs,nwatlats,Nrphys
c915144f5d Andr*0019       integer ntimesoz,ntimesq
                0020       _RL o3(nozlats,nozlevs,ntimesoz)
                0021       _RL water(nwatlats,nwatlevs,ntimesq)
328acc481c Andr*0022       _RL ozlats(nozlats), ozlevs(nozlevs)
                0023       _RL watlats(nwatlats), watlevs(nwatlevs)
                0024       _RL pressure(Nrphys),methane(Nrphys),n2o(Nrphys)
                0025       _RL co2,cfc11,cfc12,cfc22
                0026
a456aa407c Andr*0027       _RL getcon
9a6b9d7b6d Andr*0028       integer kqz,koz
                0029       call mdsfindunit( kqz, myThid )
                0030       call mdsfindunit( koz, myThid )
328acc481c Andr*0031
c915144f5d Andr*0032       call read_qz (kqz,water,watlats,watlevs,nwatlats,nwatlevs,ntimesq)
                0033       call read_oz (koz,o3,ozlats,ozlevs,nozlats,nozlevs,ntimesoz)
9a6b9d7b6d Andr*0034
328acc481c Andr*0035       call get_methane_n2o (pressure,Nrphys,n2o,methane)
9a6b9d7b6d Andr*0036
328acc481c Andr*0037       co2   = getcon('CO2'  )*1.e-6
                0038       cfc11 = getcon('CFC11')*1.e-9
                0039       cfc12 = getcon('CFC12')*1.e-9
                0040       cfc22 = getcon('CFC22')*1.e-9
                0041
                0042       RETURN
                0043       END
                0044
                0045       subroutine read_qz (ku,qz,lats,levs,nlat,nlev,ntime)
                0046 C***********************************************************************
                0047 C  PURPOSE
                0048 C     To Read Stratospheric Moisture Data
                0049 C
                0050 C  ARGUMENTS   DESCRIPTION
                0051 C     ku ...... Unit to Read  Moisture Data
                0052 C     qz ...... Stratospheric Moisture Data
                0053 C     lats .... Stratospheric Moisture Data Latitudes (degrees)
                0054 C     levs .... Stratospheric Moisture Data Levels    (mb)
                0055 C     nlat .... Number of ozone latitudes
                0056 C     nlev .... Number of ozone levels
                0057 C     ntime ... Number of ozone time values
                0058 C
                0059 C***********************************************************************
                0060
                0061       implicit none
                0062       integer  ku,nlat,nlev,ntime
                0063
9a6b9d7b6d Andr*0064       _RL qz(nlat,nlev,ntime)
a20b61c7ed Andr*0065       _RL lats(nlat)
                0066       _RL levs(nlev)
328acc481c Andr*0067
e26e993248 Andr*0068       integer time0
328acc481c Andr*0069       integer lat
                0070       integer lev
                0071
a456aa407c Andr*0072       _RL voltomas
328acc481c Andr*0073       parameter ( voltomas = 0.622e-6 )
                0074
f94be9654c Andr*0075       open(ku,file='data.sage',form='formatted')
328acc481c Andr*0076       rewind ku
                0077
                0078 c Set Moisture Data Latitudes
                0079 c ---------------------------
                0080       do   lat = 1,nlat
                0081       lats(lat) = -85. + (lat-1)*10.
                0082       enddo
                0083
                0084 c Read Moisture Pressure Levels
                0085 c -----------------------------
                0086       read(ku,1000) (levs(lev),lev=1,nlev)
                0087
                0088 c Read Moisture Amounts by Month and Level
                0089 c ----------------------------------------
e26e993248 Andr*0090       do time0=1,ntime
328acc481c Andr*0091       read (ku,1001)
                0092       do  lat=1,nlat
e26e993248 Andr*0093       read(ku,1000) (qz(lat,lev,time0),lev=1,nlev)
328acc481c Andr*0094       enddo
                0095       enddo
                0096
                0097 c Convert from Volume Mixing Ratio to Mass Mixing Ratio
                0098 c -----------------------------------------------------
e26e993248 Andr*0099       do time0 = 1,ntime
328acc481c Andr*0100       do  lev = 1,nlev
                0101       do  lat = 1,nlat
e26e993248 Andr*0102       qz(lat,lev,time0) = qz(lat,lev,time0)*voltomas
328acc481c Andr*0103       enddo
                0104       enddo
                0105       enddo
                0106
                0107  1000 format (3(5x,7(2x,f6.1)/))
                0108  1001 format (1x)
                0109       return
                0110       end
                0111
                0112       subroutine read_oz (ku,oz,lats,levs,nlat,nlev,ntime)
                0113 C***********************************************************************
                0114 C  PURPOSE
                0115 C     To Read Ozone Value
                0116 C
                0117 C  ARGUMENTS   DESCRIPTION
                0118 C     ku ...... Unit to Read Ozone Data
                0119 C     oz ...... Ozone Data
                0120 C     lats .... Ozone Data Latitudes (degrees)
                0121 C     levs .... Ozone Data Levels    (mb)
                0122 C     nlat .... Number of ozone latitudes
                0123 C     nlev .... Number of ozone levels
                0124 C     ntime ... Number of ozone time values
                0125 C
                0126 C***********************************************************************
                0127
                0128       implicit none
                0129       integer  ku,nlat,nlev,ntime
                0130
                0131       _RL   oz(nlat,nlev,ntime)
31c82eaf90 Andr*0132       real*4 o3(nlat)
a20b61c7ed Andr*0133       _RL lats(nlat)
                0134       _RL levs(nlev)
328acc481c Andr*0135
e26e993248 Andr*0136       integer time0
328acc481c Andr*0137       integer lat
                0138       integer lev
                0139       integer nrec
                0140
a20b61c7ed Andr*0141       _RL plevs(34)
328acc481c Andr*0142       data plevs/ 0.003, 0.005, 0.007, 0.01, 0.015, 0.02, 0.03, 0.05,
                0143      .            0.07, 0.1, 0.15, 0.2, 0.3, 0.5, 0.7, 1.0, 1.5, 2.0,
                0144      .            3.0, 5.0, 7.0, 10.0, 15.0, 20.0, 30.0, 50.0, 70.0,
                0145      .            100.0, 150.0, 200.0, 300.0, 500.0, 700.0, 1000.0 /
                0146
                0147 c Set Ozone Data Latitudes
                0148 c ------------------------
                0149       do   lat = 1,nlat
                0150       lats(lat) = -90. + (lat-1)*5.
                0151       enddo
                0152
                0153 c Set Ozone Data Levels
                0154 c ------------------------
                0155       do   lev = 1,nlev
                0156       levs(lev) = plevs(lev)
                0157       enddo
                0158
                0159 c Read Ozone Amounts by Month and Level
                0160 c -------------------------------------
                0161       close (ku)
f94be9654c Andr*0162       open(ku,file='data.gcmo3',form='unformatted',access='direct',
c915144f5d Andr*0163      .                                                 recl=nlat*4)
328acc481c Andr*0164
e26e993248 Andr*0165       do time0=1,ntime
328acc481c Andr*0166       do lev=1,nlev
63416ca6a5 Andr*0167 C Note: 2 quantities in Ozone Dataset
e26e993248 Andr*0168       nrec = lev+(time0-1)*nlev*2
328acc481c Andr*0169       read(ku,rec=nrec) o3
67dfde572d Andr*0170 #if defined( _BYTESWAPIO )
46f327a071 Andr*0171       call mds_byteswapr4(nlat,o3)
67dfde572d Andr*0172 #endif
328acc481c Andr*0173            do lat=1,nlat
e26e993248 Andr*0174            oz(lat,nlev-lev+1,time0) = o3(lat)
328acc481c Andr*0175            enddo
                0176       enddo
                0177       enddo
                0178
                0179       close (ku)
                0180       return
                0181       end
                0182
                0183       subroutine get_methane_n2o (pres,Nrphys,n2o,methane)
                0184 C***********************************************************************
                0185 C  PURPOSE
                0186 C     Compute methane and n2o
                0187 C
                0188 C  ARGUMENTS   DESCRIPTION
                0189 C
                0190 C***********************************************************************
                0191 C*         Climatological Annual and Global Mean Height Data           *
                0192 C***********************************************************************
                0193
                0194       implicit none
                0195
                0196       integer Nrphys
9a6b9d7b6d Andr*0197       _RL n2o(Nrphys),methane(Nrphys)
328acc481c Andr*0198       _RL pres(Nrphys)
a20b61c7ed Andr*0199       _RL hght(Nrphys), slope,pr1,pr2,hpr1,hpr2
328acc481c Andr*0200       integer L,L1,L2,lup,ldn
                0201
a20b61c7ed Andr*0202       _RL plevc (46), plevz(46)
                0203       _RL hghtc (46), hghtz(46)
328acc481c Andr*0204
                0205       data plevc /1000.00, 975.00, 950.00, 925.00, 900.00,
                0206      .             875.00, 850.00, 825.00, 800.00, 750.00,
                0207      .             700.00, 650.00, 600.00, 550.00, 500.00,
                0208      .             450.00, 400.00, 350.00, 300.00, 250.00,
                0209      .             200.00, 150.00, 100.00,  70.00,  50.00,
                0210      .              40.00,  30.00,  20.00,  10.00,   7.00,
                0211      .               5.00,   4.00,   3.00,   2.00,   1.00,
                0212      .               0.70,   0.50,   0.40,   0.30,   0.20,
                0213      .               0.10,   0.07,   0.05,   0.04,   0.03,
                0214      .               0.02 /
                0215
                0216       data hghtc/ 0.128733 , 0.316985 , 0.528275 , 0.749515 , 0.976471 ,
                0217      .            1.208910 , 1.446800 , 1.690980 , 1.941630 , 2.463530 ,
                0218      .            3.016200 , 3.603490 , 4.229410 , 4.899870 , 5.622320 ,
                0219      .            6.405940 , 7.263450 , 8.211920 , 9.275540 , 10.49150 ,
                0220      .            11.92420 , 13.70200 , 16.12980 , 18.24120 , 20.26480 ,
                0221      .            21.63100 , 23.41250 , 25.96570 , 30.45890 , 32.85240 ,
                0222      .            35.17360 , 36.75040 , 38.82900 , 41.84600 , 47.15580 ,
                0223      .            49.90100 , 52.46230 , 54.13890 , 56.26340 , 59.17640 ,
                0224      .            63.89980 , 66.20240 , 68.29210 , 69.63550 , 71.32330 ,
                0225      .            73.62110 /
                0226
                0227       do L=1,46
                0228       plevz(L) = plevc(47-L)
                0229       hghtz(L) = hghtc(47-L)
                0230       enddo
                0231
                0232 C **********************************************************************
                0233 C                     Interpolate Heights to Model Pressures        ****
                0234 C **********************************************************************
                0235
                0236       do L2 = 1,Nrphys
                0237
                0238          do L1 = 1,46
                0239             if( plevz(L1).gt.pres(L2) ) then
                0240                 if( L1.eq.1 ) then
                0241                     lup = 1
                0242                     ldn = 2
                0243                 else
                0244                     lup = L1-1
                0245                     ldn = L1
                0246                 endif
                0247                 goto 10
                0248             endif
                0249          enddo
                0250          lup = 45
                0251          ldn = 46
                0252
                0253    10 continue
                0254        pr1 = plevz(lup)
                0255        pr2 = plevz(ldn)
                0256       hpr1 = hghtz(lup)
                0257       hpr2 = hghtz(ldn)
                0258
                0259       slope = ( hpr1-hpr2 )/( pr1-pr2 )
                0260       hght(L2) = hpr2 + ( pres(L2)-pr2 )*slope
                0261
                0262       enddo
                0263
                0264 C **********************************************************************
                0265 C  Set the profiles of N2O and CH4 based on Bresser and Pawson 1996 ****
                0266 C **********************************************************************
                0267
                0268       do L = 1,Nrphys
                0269        if( hght(L).gt.26. ) then
                0270         n2o(L) = 120.* exp( (26.- hght(L)) / 5.69 ) * 1.e-9
                0271        else if( hght(L).gt.16. ) then
                0272         n2o(L) = 307.* exp( (16.- hght(L)) /10.47 ) * 1.e-9
                0273        else
                0274         n2o(L) = 307.e-9
                0275        endif
                0276       enddo
                0277
                0278       do L = 1,Nrphys
                0279        if( hght(L).gt.55. ) then
                0280         methane(L) = 0.2e-6
                0281        else if( hght(L).gt.14. ) then
                0282         methane(L) = 1.7* exp( (14.- hght(L)) /19.16 ) * 1.e-6
                0283        else
                0284         methane(L) = 1.7e-6
                0285        endif
                0286       enddo
                0287
                0288       return
                0289       end