C**********************************************************                  MORC.2     
C             SUBROUTINE RADN                                                MORC.3     
C**********************************************************                  MORC.4     
      SUBROUTINE RADN(TROPHT)                                                MORC.5     
C                                                                            MORC.6     
C     RADIATION SCHEME. MORCRETTE PMF 13/5/97                                MORC.7     
C     It passes the pressure of the                                          MORC.8     
C     full sigma levels and the surface                                      MORC.9     
C     to the Radiation scheme                                                MORC.10    
C     (So the morcrette code has one more level than NL)                     MORC.11    
C     water vapour values come from QG and QSTAR                             MORC.12    
C     and temperatures from TG and TSTAR                                     MORC.13    
C                                                                            MORC.14    
C     Dingmin's changes to morcrette                                         MORC.15    
C     combine bottom levels                                                  MORC.16    
C     and set bottom level heating rate to zero as                           MORC.17    
C     this is now the ground and is not a layer                              MORC.18    
C                                                                            MORC.19    
C                                                                            PARAM1.2     
C     Determines model resolution                                            PARAM1.3     
C                                                                            PARAM1.4     
      PARAMETER(NN=21,MM=21,NHEM=2,NL=5,MOCT=1,MG=64,JG=16,NWJ2=121          NETCDF.35    
     P         ,NCRAY=64,JGL=JG,NTRAC=1,NLEVRF=1)                            NETCDF.36    
                                                                             PARAM1.7     
C                                                                            PARAM1.8     
C                                                                            PARAM2.2     
C     Sets basic constants, especially those needed for array dimensions     PARAM2.3     
C                                                                            PARAM2.4     
      PARAMETER(MH=2,PI=3.14159265359,PI2=2.0*PI                             PARAM2.5     
     +,NNP=NN+1,MGPP=MG+2,JGP=JG+1,JGG=JG*NHEM,JGGP=JGG+1,MJP=NWJ2+NWJ2      PARAM2.6     
     +,NLM=NL-1,NLP=NL+1,NLPP=NL+2,NLA=NL+3,NLB=NL+4,NL2=NL*NL               PARAM2.7     
     +,IDA=(MG+MG+MG)/2+1,IDB=NWJ2*NL,IDC=IDB+IDB,IDD=MGPP*NL                PARAM2.8     
     +,IDE=NL2*NN,IDF=NCRAY*(MG+1),IDG=JG*NL,IDH=JG*MG                       PARAM2.9     
     +,IDI=NNP/2,IDJ=IDI*IDI,IDK=NL*IDI,IDL=MGPP/2,IDM=NNP/2,IDN=IDM*NL      PARAM2.10    
     +,NWW=1+(MM-1)/MOCT)                                                    PARAM2.11    
      PARAMETER(IGA=NWJ2*NHEM,IGB=IDB*NHEM,IGC=MGPP*NHEM,IGD=IDD*NHEM        PARAM2.12    
     +,IGG=IDG*NHEM,IGL=IDL*NHEM,IGM=IDM*NHEM,IGN=IDN*NHEM                   PARAM2.13    
     +,IGO=IGA+IGA,IGP=IGB+IGB,NFTWG=(5+NTRAC)*NL+3                          PARAM2.14    
     +,NFTGW=(6+3*NTRAC)*NL+2,NFTGD=(3+NTRAC)*NL,NLTR=NL*NTRAC)              PARAM2.15    
C     Number of 2D (surface) output fields. This value is                    PARAM2.16    
C     Doubled due to averaged and instantaneous fields.                      PARAM2.17    
      PARAMETER (N2DFLD=21,NGRPAD=N2DFLD*2*IGC)                              PARAM2.18    
C                                                                            PARAM2.19    
C                                                                            BLANK.2     
C     Basic planetary parameters for run plus information about              BLANK.3     
C     vertical grid structure                                                BLANK.4     
C                                                                            BLANK.5     
C     Note that RD and GASCON are identical and CPD is set from RD,AKAP.     BLANK.6     
      COMMON        SQ(NNP),RSQ(NNP),SIGMAH(NLM),SIGMA(NL)                   BLANK.7     
     +              ,T01S2(NLM),T0(NL),ALPHA(NL),DSIGMA(NL),RDSIG(NL)        BLANK.8     
     +              ,TKP(NL),C(NL2),SQH(NNP)                                 BLANK.9     
     +              ,MF,MFP,JZF,NF,NFP                                       BLANK.10    
     +              ,AKAP,GA,GASCON,RADEA,WW,P0,PFAC,EZ,AIOCT                BLANK.11    
     +              ,RD,RV,CPD,CLATNT                                        BLANK.12    
     +              ,LRSTRT,LSHORT,LTVEC,LSTRETCH                            BLANK.13    
     +              ,LFLUX                                                   BLANK.14    
     +              ,LBALAN,LRESTIJ                                          BLANK.15    
     +              ,LCLIM, LPERPET, L22L,LOROG ,LCSFCT                      BLANK.16    
     +              ,LNOISE                                                  BLANK.17    
      COMPLEX EZ,AIOCT                                                       BLANK.18    
      LOGICAL LRSTRT,LSHORT,LTVEC,LSTRETCH,LBALAN,LRESTIJ                    BLANK.19    
     +       ,LFLUX,LNOISE                                                   BLANK.20    
     +       ,LCLIM, LPERPET, L22L,LOROG,LCSFCT                              BLANK.21    
C                                                                            BLANK.22    
C                                                                            LEGAU.2     
C     Legendre polynomials and information about gaussian latitudes          LEGAU.3     
C                                                                            LEGAU.4     
      COMMON/LEGAU/ ALPJ(MJP),DALPJ(MJP)                                     LEGAU.5     
     +              ,ALP(NWJ2,2,JGL),DALP(NWJ2,2,JGL)                        LEGAU.6     
     +              ,RLP(NWJ2,2,JGL),RDLP(NWJ2,2,JGL)                        LEGAU.7     
     +              ,SI(JGG),CS(JGG),SISQ(JGG),CSSQ(JGG),SECSQ(JGG)          LEGAU.8     
     +              ,ALAT(JGG),GWT(JGG),AW(JGG),JH,JL,JINC                   LEGAU.9     
C                                                                            LEGAU.10    
C                                                                            GRIDPP.2     
C     Array ordering in GRIDP must correspond to that in SPECTR.             GRIDPP.3     
C     Real arrays: multi-level arrays are 2-dimensional.                     GRIDPP.4     
C     the variables have been renamed to coincide with                       GRIDPP.5     
C     variable names in bgcm5 DGRMLT                                         GRIDPP.6     
C                                                                            GRIDPP.7     
C swapped round VLNG and UNLG to check                                       GRIDPP.8     
      COMMON/GRIDP/ CHIG(IGC,NL),SFG(IGC,NL),UG(IGC,NL),VG(IGC,NL)           GRIDPP.9     
     :              ,TTVD(IGC,NL),QTVD(IGC,NL),TG(IGC,NL)                    GRIDPP.10    
     :              ,TRAG(IGC,NL,NTRAC)                                      GRIDPP.11    
     :              ,PLG(IGC),TYBL(IGC),TXBL(IGC)                            GRIDPP.12    
     :              ,SPG(IGC),VPG(IGC),TTRD(IGC,NL)                          GRIDPP.13    
     :              ,TNLG(IGC,NL),TRANLG(IGC,NL,NTRAC),UNLG(IGC,NL)          GRIDPP.14    
     :              ,VNLG(IGC,NL),TTLR(IGC,NL),UTRAG(IGC,NL,NTRAC)           GRIDPP.15    
     :              ,TTCR(IGC,NL),VTRAG(IGC,NL,NTRAC)                        GRIDPP.16    
     :              ,UTVD(IGC,NL),VTVD(IGC,NL)                               GRIDPP.17    
     :         ,ASSBL(IGC),ASHBL(IGC),ASLBL(IGC),ARRCR(IGC),ARRLR(IGC)       GRIDPP.18    
     :         ,arflux(igc,6),asfld(igc,6),acld(igc,4)                       GRIDPP.19    
     :         ,SSBL(IGC),SHBL(IGC),SLBL(IGC),RRCR(IGC),RRLR(IGC)            GRIDPP.20    
     :         ,rflux(igc,6),sfld(igc,6),cld(igc,4)                          GRIDPP.21    
C                                                                            GRIDPP.22    
C                                                                            BATS.2     
C     Constant arrays and variables associated with time and vertical        BATS.3     
C     differencing. Also counters.                                           BATS.4     
C                                                                            BATS.5     
      COMMON/BATS/  BM1(IDE),AK(NNP),AQ(NL2),G(NL2),TAU(NL2)                 BATS.6     
     +              ,KOUNT,KITS,KSTART,KTOTAL,KRUN,BEGDAY,ITSPD              BATS.7     
     +              ,DELT,DELT2,CV,CG,CT,CQ,PNU,PNU2,PNU21                   BATS.8     
     +              ,NTRACO,CTRA(NTRAC),KOLOUR(NTRAC),RGG(NL2)               BATS.9     
     +              ,BEGDOY,DOY                                              BATS.10    
C                                                                            BATS.11    
      COMMON/PHYS/  LBL,LVD,LCR,LLR,LRD,LCUBM,LCBADJ                         PHYSCA.2     
     :              ,TSTAR(IGC,JG),QSTAR(IGC,JG),FRAD(JG,NHEM)               PHYSCA.3     
     :              ,TSTARO(IGC,JG),TDEEPO(IGC,JG),smstar(igc,jg)            PHYSCA.4     
     :              ,tdeep(igc,jg),hsnow(igc,jg),sqstar(igc,jg)              PHYSCA.5     
     :              ,SALB(IGC,JG),SBAL(IGC,JG),BLCD(IGC)                     PHYSCA.6     
     :              ,SVEGE(IGC,JG),CD,DRAG,BLVAD,BLA,BLRH,BLVB(IGC)          PHYSCA.7     
     :              ,AKVV,AKTV,AKQV,ESCONA,ESCONB,EPSIQ,CTQ,CCC              PHYSCA.8     
     : ,ctqi,sdsn,shcs,shcsp,shcsn,skse,sksn,slhf,sd1,sd2,sdw                PHYSCA.9     
     :        ,ssmc,sdsnd,LSL,sasnow,saice,shsstar,shsmax                    PHYSCA.10    
     :     ,LOC,SHCO,SHCI,LNOICE,LOLDBL,LCOND,LNNSK,ITSLL,ITSLO              PHYSCA.11    
     :              ,CCR,RCON,DTBUOY,TSLA,TSLB,TSLC,TSLD,CUT1,CUT2           PHYSCA.12    
     :              ,NLCR,NCUTOP,CURHM,AKTC,AKQC,CUBMT,CBADJT,CBADJP         PHYSCA.13    
     :              ,SKAP(NL),SK(NLM),FWS(NL),CLR(NL),FB(NLM)                PHYSCA.14    
     :              ,TTDC(NL),QTDC(NL),TTMC(NL),QTMC(NL),TC(NL),QC(NL)       PHYSCA.15    
     :              ,CTCR(NL,NHEM),CTLR(NL,NHEM),NAVRD,NAVWT,DELT2C          PHYSCA.16    
      LOGICAL LBL,LVD,LCR,LLR,LRD,LCUBM,LCBADJ,LSL,LOC                       PHYSCA.17    
     :       ,LNOICE,LOLDBL,LCOND,LNNSK                                      PHYSCA.18    
C                                                                            PHYSCA.19    
C                                                                            XTRAP.2     
C     Setup moisture variables by equivilencing them to                      XTRAP.3     
C     Tracer No. 1                                                           XTRAP.4     
C                                                                            XTRAP.5     
      REAL QG(IGC,NL),QNLG(IGC,NL),QTLR(IGC,NL),QTCR(IGC,NL)                 XTRAP.6     
      EQUIVALENCE (QG(1,1),TRAG(1,1,1)) , (QNLG(1,1),TRANLG(1,1,1))          XTRAP.7     
     *            ,(QTLR(1,1),UTRAG(1,1,1)),(QTCR(1,1),VTRAG(1,1,1))         XTRAP.8     
C                                                                            XTRAP.9     
       COMMON/CPIERS/ICFLAG(IGC,5,2),CFRAC(IGC,5),PNET(IGC,JG)               CPIERS.2     
     :               ,SNET(IGC,JG),RRFLUX(IGC,JG,6)                          CPIERS.3     
     :               ,TTSW(IGC,NL),TTLW(IGC,NL)                              CPIERS.4     
       COMMON/GSG/GSG(IGC,JG)                                                OROG_G.2     
       REAL htnet                                                            RADHT.2     
       COMMON /RADHT/ HTNET(NHEM,JG,MG,NL)                                   RADHT.3     
       REAL TAVE(IGP)                                                        RADHT.4     
c                                                                            MORC.31    
c Climatology ozone and water values. '1' is current month, '2' is           MORC.32    
c the next month. (SMR 30-05-97)                                             MORC.33    
c 11-03-97 i.e. nl=15 read in from climatology. Interpolated to              MORC.34    
c nl.gt.15 in subroutine interpo3.                                           MORC.35    
                                                                             MORC.36    
      real o3clim1(15,mg,nhem,jg),h2oclim1(15,mg,nhem,jg)                    MORC.37    
      real o3clim2(15,mg,nhem,jg),h2oclim2(15,mg,nhem,jg)                    MORC.38    
                                                                             MORC.39    
                                                                             MORC.40    
                                                                             MORC.41    
c Morcrette profiles (have extra level).                                     MORC.42    
                                                                             MORC.43    
      REAL PR(NL+1),T(NL+1),h2o(nl+1),o3(nl+1),htlw(nl+1),htsw(nl+1)         MORC.44    
      real adum                                                              MORC.45    
                                                                             MORC.46    
      real o3mod1(nl,mg,nhem,jg)    ! ozone interpolated to model            MORC.47    
                                    ! levels. Current month.                 MORC.48    
      real h2omod1(nl,mg,nhem,jg)   ! water interpolated to model            MORC.49    
                                    ! levels. Current month.                 MORC.50    
      real o3mod2(nl,mg,nhem,jg)    ! ozone interpolated to model            MORC.51    
                                    ! levels. Next month.                    MORC.52    
      real h2omod2(nl,mg,nhem,jg)   ! water interpolated to model            MORC.53    
                                    ! levels. Next month.                    MORC.54    
      real tropht(mg,nhem,jg)       ! tropopause height.                     MORC.55    
                                                                             MORC.56    
      real rdummy(15)               ! dummy profile for skipping to          MORC.57    
                                    ! the right part of the                  MORC.58    
                                    ! climatology                            MORC.59    
                                                                             MORC.60    
                                                                             MORC.61    
      integer ifirst                ! If =1, first time reading o3           MORC.62    
                                    ! and h2o (2 months' worth).             MORC.63    
                                                                             MORC.64    
      character*31 o3file1, h2ofile1, o3file2, h2ofile2                      MORC.65    
                                    ! Filenames of o3 and h2o                MORC.66    
                                    ! climatology read in. 1=current,        MORC.67    
                                    ! 2=next.                                MORC.68    
      integer cmth                  ! Current month counter. = month         MORC.69    
                                    ! number +1 'cos have 13 months          MORC.70    
                                    ! (dec in twice)                         MORC.71    
      character*3 mn(13)            ! Month names. 13 to enable              MORC.72    
                                    ! wrap-around.                           MORC.73    
      integer mnlen(13)             ! Month lengths in days.                 MORC.74    
      real amfrac                   ! fraction through month                 MORC.75    
      integer ichange               ! =1 when in process of month change     MORC.76    
*---------------------------                                                 MORC.77    
      integer ifirstcol             ! =1 first time through column           MORC.78    
                                    ! calculation (open new file).           MORC.79    
      real p0                                                                MORC.80    
      real ps                       ! sfc pressure (used in                  MORC.81    
                                    ! interpolation from climatology         MORC.82    
                                    ! to model).                             MORC.83    
      integer im                    ! Pointer for array plg (for             MORC.84    
                                    ! getting sfc pressure).                 MORC.85    
                                                                             MORC.86    
      integer ic(4,2)    !cloud positions (deep,h,m,l;bottom,top)            MORC.87    
      real cf(4,2)       !cloud fraction,lwp deep,high,mid,low               MORC.88    
C     Array to hold fluxes at top and bottom of atmosphere                   MORC.89    
C     1st index - flux 1=SW, 2=LW                                            MORC.90    
C     2nd index - Direction 1=DN, 2=UP                                       MORC.91    
C     3rd index - Where 1=TOP, 2=SURFACE                                     MORC.92    
      real fluxes(2,2,2)                                                     MORC.93    
                                                                             MORC.94    
      save                          ! Want to keep things like dcompl.       MORC.95    
                                                                             MORC.96    
      data mn/'dec','jan','feb','mar','apr','may','jun','jul',               MORC.97    
     &          'aug','sep','oct','nov','dec'/                               MORC.98    
                                                                             MORC.99    
                                                                             MORC.100   
      DATA IFIRST/1/                                                         MORC.101   
      data ifirstcol/1/                                                      MORC.102   
                                                                             MORC.103   
      RHSCL=288.0*GASCON/GA                                                  MORC.104   
      CHRF=86400.*WW*CT  ! factor to non-dimensionalise heating rates        MORC.105   
c  Skipping part set up here.                                                MORC.106   
c  Radiation scheme only called every nskip longitudes                       MORC.107   
c  nskip must divide exactly into mg for longitude.                          MORC.108   
c  ntstep is the number of timesteps to skip.                                MORC.109   
                                                                             MORC.110   
      ntstep=itspd       ! N.B. No point any more frequent as                MORC.111   
                         ! Morcrette doing diurnal averages anyway.          MORC.112   
      nskip=4                                                                MORC.113   
                                                                             MORC.114   
      if (jg.eq.16.and.mg.eq.64) then   ! Crude variation of nskip           MORC.115   
                                        ! with latitude.                     MORC.116   
        if (jh.eq.1) nskip=32                                                MORC.117   
        if (jh.eq.2) nskip=16                                                MORC.118   
        if (jh.eq.3.or.jh.eq.4.or.jh.eq.5) nskip=8                           MORC.119   
        if (jh.eq.13.or.jh.eq.14.or.jh.eq.15.or.jh.eq.16) nskip=2            MORC.120   
      endif                                                                  MORC.121   
                                                                             MORC.122   
c                                                                            MORC.123   
C this sets the non-dimensional gridpoint temperature tendency (TTRD)        MORC.124   
C to get to TTRD from K/day divide the K/day heating rate by                 MORC.125   
C     (86400*WW*CT)                                                          MORC.126   
C     TTRD(non-dim)=HTRT(k/day)/(86400*WW*CT)                                MORC.127   
C                                                                            MORC.128   
c                                                                            MORC.129   
C Reads in bogus ozone and water vapour on first call                        MORC.130   
C other data such as albedo and number of levels is in the include file      MORC.131   
C                                                                            MORC.132   
                                                                             MORC.133   
                                                                             MORC.134   
c ---------- Read in required ozone and water vapour from climatology        MORC.135   
                                                                             MORC.136   
                                                                             MORC.137   
      IF (ifirst.eq.1) then           ! First time through (read in          MORC.138   
                                      ! two months of o3 and h2o if          MORC.139   
                                      ! lperpet .false. or just one if       MORC.140   
                                      ! lperpet .true.).                     MORC.141   
                                                                             MORC.142   
       IF (JH.EQ.1) THEN                                                     MORC.143   
        CALL CALNDR(DOY,MTH1,AMFRAC)                                         MORC.144   
        CMTH=MTH1+1                                                          MORC.145   
                                      ! in first.                            MORC.146   
        IF(cmth.eq.13) cmth=1        ! Wrap-around to January.               MORC.147   
                                                                             MORC.148   
                                                                             MORC.149   
c ------------------------- Make the filenames and open the files:           MORC.150   
c                           N.B. Read in from global climatology.            MORC.151   
                                                                             MORC.152   
        write(o3file1,'(a,i3.3,a,i3.3,a)')                                   MORC.153   
     &  'ozone_'//mn(cmth)//'_jg',jg,'_mg',mg,'_nhem2.dat'                   MORC.154   
        write(h2ofile1,'(a,i3.3,a,i3.3,a)')                                  MORC.155   
     &  'water_'//mn(cmth)//'_jg',jg,'_mg',mg,'_nhem2.dat'                   MORC.156   
                                                                             MORC.157   
        write(o3file2,'(a,i3.3,a,i3.3,a)')                                   MORC.158   
     &  'ozone_'//mn(cmth+1)//'_jg',jg,'_mg',mg,'_nhem2.dat'                 MORC.159   
        write(h2ofile2,'(a,i3.3,a,i3.3,a)')                                  MORC.160   
     &  'water_'//mn(cmth+1)//'_jg',jg,'_mg',mg,'_nhem2.dat'                 MORC.161   
                                                                             MORC.162   
        open (47,file='climdata/'//o3file1,status='old')                     MORC.163   
        open (49,file='climdata/'//h2ofile1,status='old')                    MORC.164   
        open (51,file='climdata/'//o3file2,status='old')                     MORC.165   
        open (53,file='climdata/'//h2ofile2,status='old')                    MORC.166   
       ENDIF                                                                 MORC.167   
                                                                             MORC.168   
       j =jh                                                                 MORC.169   
       DO 796 ihem=1,nhem      ! start of loop over hemispheres              MORC.170   
        DO 798 i=1,mg       ! start of loop over longitude                   MORC.171   
         read(47,*)(o3clim1(l,i,ihem,j),l=1,15)                              MORC.172   
         read(49,*)(h2oclim1(l,i,ihem,j),l=1,15)                             MORC.173   
         read(51,*)(o3clim2(l,i,ihem,j),l=1,15)                              MORC.174   
         read(53,*)(h2oclim2(l,i,ihem,j),l=1,15)                             MORC.175   
                                                                             MORC.176   
c ----------------------------- Interpolate to model vertical grid.          MORC.177   
                                                                             MORC.178   
         im=i+(ihem-1)*mgpp                                                  MORC.179   
         ps=1.0E5                                                            MORC.180   
                                                                             MORC.181   
         call interpf(h2oclim1(1,i,ihem,j),h2omod1(1,i,ihem,j),ps)           MORC.182   
         call O3INTERP(o3clim1(1,i,ihem,j),o3mod1(1,i,ihem,j),ps)            MORC.183   
         call interpf(h2oclim2(1,i,ihem,j),h2omod2(1,i,ihem,j),ps)           MORC.184   
         call O3INTERP(o3clim2(1,i,ihem,j),o3mod2(1,i,ihem,j),ps)            MORC.185   
                                                                             MORC.186   
  798   CONTINUE            ! end of loop over longitude                     MORC.187   
                                                                             MORC.188   
        IF (nhem.eq.1) then     ! skip the other hemisphere of the           MORC.189   
                                 ! climatology (which is global).            MORC.190   
         DO i=1,mg                                                           MORC.191   
          read(47,*)(rdummy(l),l=1,15)                                       MORC.192   
          read(49,*)(rdummy(l),l=1,15)                                       MORC.193   
          read(51,*)(rdummy(l),l=1,15)                                       MORC.194   
          read(53,*)(rdummy(l),l=1,15)                                       MORC.195   
         ENDDO                                                               MORC.196   
        ENDIF                                                                MORC.197   
                                                                             MORC.198   
  796  CONTINUE              ! end of loop over hemispheres                  MORC.199   
                                                                             MORC.200   
                                                                             MORC.201   
c ------------- Close the files:                                             MORC.202   
       IF (JH.EQ.JG) THEN                                                    MORC.203   
        IFIRST=0                                                             MORC.204   
                                                                             MORC.205   
        close(47)                                                            MORC.206   
        close(49)                                                            MORC.207   
        close(51)                                                            MORC.208   
        close(53)                                                            MORC.209   
!  Open file for column amount calculations.                                 MORC.210   
        IF (ifirstcol.eq.1) then                                             MORC.211   
         ifirstcol=0                                                         MORC.212   
         open(unit=55,file='column.dat',status='new')                        MORC.213   
        ENDIF                                                                MORC.214   
                                                                             MORC.215   
        write(55,*)'******** In ifirst loop'                                 MORC.216   
        write(55,*)'******** Start month is: ',mn(cmth)                      MORC.217   
        write(55,*)'**********************************'                      MORC.218   
        call colamt(o3mod1,h2omod1)                                          MORC.219   
       ENDIF                                                                 MORC.220   
                                                                             MORC.221   
      ENDIF       ! end of first time bit (reading in and                    MORC.222   
                  ! interpolating o3 and h2o).                               MORC.223   
                                                                             MORC.224   
c -------------------- If it's time to change month of o3 (and h2o)          MORC.225   
c                      for a time-varying o3 run.                            MORC.226   
c                                                                            MORC.227   
      IF (JH.EQ.1) CALL CALNDR(DOY,MTH2,AMFRAC)                              MORC.228   
      IF (.not. lperpet) then                                                MORC.229   
       IF (MTH1.NE.MTH2) THEN                                                MORC.230   
c  begin smrdebug                                                            MORC.231   
        IF (JH.EQ.1) THEN                                                    MORC.232   
         ifirstdebug=1                                                       MORC.233   
                                                                             MORC.234   
         write(55,*)'**************************************'                 MORC.235   
         write(55,*)'Start of ''time to change ozone'' block: '              MORC.236   
         write(55,*)'kount: ',kount                                          MORC.237   
         write(55,*)'kount/itspd i.e. no. of days: ',kount/itspd             MORC.238   
         write(55,*)'Month changing from ',mn(cmth),' to ',                  MORC.239   
     &               mn(cmth+1)                                              MORC.240   
         write(55,*)'*********************************'                      MORC.241   
c end smrdebug                                                               MORC.242   
                                                                             MORC.243   
         CMTH=MTH2+1                                                         MORC.244   
                                                                             MORC.245   
         IF (cmth.eq.13) cmth=1      ! Wrap-around to January.               MORC.246   
                                                                             MORC.247   
! Copy 'next' to 'current'.                                                  MORC.248   
         do j=1,jg         ! loop over latitude                              MORC.249   
          do ihem=1,nhem     !           hemispheres                         MORC.250   
           do i=1,mg       !           longitude                             MORC.251   
            do l=1,nl      !           height                                MORC.252   
             o3mod1(l,i,ihem,j)=o3mod2(l,i,ihem,j)                           MORC.253   
             h2omod1(l,i,ihem,j)=h2omod2(l,i,ihem,j)                         MORC.254   
            enddo                                                            MORC.255   
           enddo                                                             MORC.256   
          enddo                                                              MORC.257   
         enddo                                                               MORC.258   
                                                                             MORC.259   
! And fetch new 'next':                                                      MORC.260   
         write(o3file2,'(a,i3.3,a,i3.3,a)')                                  MORC.261   
     &  'ozone_'//mn(cmth+1)//'_jg',jg,'_mg',mg,'_nhem2.dat'                 MORC.262   
         write(h2ofile2,'(a,i3.3,a,i3.3,a)')                                 MORC.263   
     &  'water_'//mn(cmth+1)//'_jg',jg,'_mg',mg,'_nhem2.dat'                 MORC.264   
                                                                             MORC.265   
         open (51,file='climdata/'//o3file2,status='old')                    MORC.266   
         open (53,file='climdata/'//h2ofile2,status='old')                   MORC.267   
        ENDIF                                                                MORC.268   
                                                                             MORC.269   
        j=jh                                                                 MORC.270   
        do ihem=1,nhem           !           hemispheres                     MORC.271   
         do i=1,mg             !           longitude                         MORC.272   
                                                                             MORC.273   
          read(51,*)(o3clim2(l,i,ihem,j),l=1,15)                             MORC.274   
          read(53,*)(h2oclim2(l,i,ihem,j),l=1,15)                            MORC.275   
                                                                             MORC.276   
          im=i+(ihem-1)*mgpp                                                 MORC.277   
          ps=1.0E5                                                           MORC.278   
                                                                             MORC.279   
          call interpf(h2oclim2(1,i,ihem,j),h2omod2(1,i,ihem,j),ps)          MORC.280   
          call O3INTERP(o3clim2(1,i,ihem,j),o3mod2(1,i,ihem,j),ps)           MORC.281   
                                                                             MORC.282   
         ENDDO                  ! end of loop over longitude                 MORC.283   
                                                                             MORC.284   
         IF (nhem.eq.1) then    ! skip the other hemisphere of               MORC.285   
                                ! the climatology                            MORC.286   
          DO i=1,mg                                                          MORC.287   
           read(51,*)(rdummy(l),l=1,15)                                      MORC.288   
           read(53,*)(rdummy(l),l=1,15)                                      MORC.289   
          ENDDO                                                              MORC.290   
         ENDIF                                                               MORC.291   
                                                                             MORC.292   
        ENDDO                  ! end of loop over hemispheres                MORC.293   
                                                                             MORC.294   
        IF (JH.EQ.JG) THEN                                                   MORC.295   
         MTH1=MTH2   ! stops loop repeating                                  MORC.296   
         close(51)                                                           MORC.297   
         close(53)                                                           MORC.298   
                                                                             MORC.299   
         call colamt(o3mod1,h2omod1)                                         MORC.300   
        ENDIF                                                                MORC.301   
       ENDIF                                                                 MORC.302   
      ENDIF                                                                  MORC.303   
c --------------------------------- Now start the radiation bit.             MORC.304   
c                                                                            MORC.305   
C loop over hemispheres                                                      MORC.306   
      IOFM=0                                                                 MORC.307   
C                                                                            MORC.308   
      DO 800 ihem=1,nhem                                                     MORC.309   
                                                                             MORC.310   
c calculates heating rates every ntstep time steps                           MORC.311   
c                                                                            MORC.312   
C                                                                            MORC.313   
       IF (mod(kount,ntstep).eq.1) then                                      MORC.314   
C                                                                            MORC.315   
C  Does do Radn scheme                                                       MORC.316   
C                                                                            MORC.317   
C                                                                            MORC.318   
C loop over longitudes for radn calculation                                  MORC.319   
      ilast=0                                                                NETCDF.51    
        DO i=1,mg                                                            MORC.320   
         im=i+iofm                                                           MORC.321   
         idocalc=0                                                           MORC.322   
         IF ((i.eq.1).or.(i-ilast.ge.nskip)) then                            MORC.323   
          idocalc=1                                                          MORC.324   
         ELSE                                                                MORC.325   
          IF (LNNSK) THEN                                                    MORC.326   
           imp=im+1                                                          MORC.327   
           IF (imp.gt.(mg+iofm)) imp=1+iofm                                  MORC.328   
           imm=im-1                                                          MORC.329   
           IF (((gsg(im,jh).gt.0.).and.(gsg(imp,jh).eq.0.)).or.              MORC.330   
     $      ((gsg(im,jh).eq.0.).and.(gsg(imp,jh).gt.0.)).or.                 MORC.331   
     $      ((gsg(im,jh).gt.0.).and.(gsg(imm,jh).eq.0.)).or.                 MORC.332   
     $      ((gsg(im,jh).eq.0.).and.(gsg(imm,jh).gt.0.))) THEN               MORC.333   
            idocalc=1                                                        MORC.334   
           ENDIF                                                             MORC.335   
          ENDIF                                                              MORC.336   
         ENDIF                                                               MORC.337   
         IF (idocalc.eq.1) then                                              MORC.338   
                                                                             MORC.339   
c ------------------------------ First set SURFACE VALUES                    MORC.340   
c Pressure, units of Pa                                                      MORC.341   
          PR(1)=PLG(im)*P0                                                   MORC.342   
c T, units of K                                                              MORC.343   
          T(1)=TSTAR(IM,JH)*CT                                               MORC.344   
c water, units of mmr                                                        MORC.345   
                                                                             MORC.346   
          IF (lclim) then             ! If using h2o from climatology.       MORC.347   
           h2o(1)=AMFRAC*h2omod2(NL,i,ihem,JH)+                              MORC.348   
     +            (1.0-AMFRAC)*h2omod1(NL,i,ihem,JH)                         MORC.349   
          ELSE                           ! using h2o from model.             MORC.350   
           h2o(1)=QSTAR(IM,JH)                                               MORC.351   
          ENDIF                                                              MORC.352   
                                                                             MORC.353   
          o3(1)=AMFRAC*o3mod2(NL,i,ihem,JH)+                                 MORC.354   
     +            (1.0-AMFRAC)*o3mod1(NL,i,ihem,JH)                          MORC.355   
                                                                             MORC.356   
c --------------------------------------- Now set rest of column.            MORC.357   
                                                                             MORC.358   
          DO LD=1,NL                     ! Start of loop over column.        MORC.359   
           L=NL+2-LD                     ! Reverse index (Morc goes          MORC.360   
                                        ! bottom up).                        MORC.361   
           PR(L)=SIGMA(LD)*PR(1)         ! Pressure                          MORC.362   
           T(L)=TG(im,ld)*CT             ! Temperature                       MORC.363   
                                                                             MORC.364   
C Water vapour, including dealing with minimum values.                       MORC.365   
                                                                             MORC.366   
           IF (lclim) then            ! Use h2o from climatology.            MORC.367   
                                                                             MORC.368   
            h2o(l)=AMFRAC*h2omod2(ld,i,ihem,JH)+                             MORC.369   
     +            (1.0-AMFRAC)*h2omod1(ld,i,ihem,JH)                         MORC.370   
                                                                             MORC.371   
            h2o(L)=MAX(6.0E-6,h2o(L))                                        MORC.372   
C  fix water vapour for levels above climatological tropopause to be         MORC.373   
C  6e-6 ppmv.                                                                MORC.374   
            IF(PR(L).LE.TROPHT(i,ihem,JH))THEN                               MORC.375   
C  remove following line if no fix required                                  MORC.376   
               h2o(L)=6.0E-6                                                 MORC.377   
            ENDIF ! end of strat. water vapour fix                           MORC.378   
                                                                             MORC.379   
           ELSE                       ! h2o from model                       MORC.380   
                                                                             MORC.381   
            h2o(L)=MAX(6.0E-6,QG(im,LD))                                     MORC.382   
                                                                             MORC.383   
C  fix water vapour for levels above climatological tropopause               MORC.384   
            IF(PR(L).LE.TROPHT(i,ihem,JH))THEN                               MORC.385   
               h2o(L)=6.0E-6                                                 MORC.386   
            ENDIF ! end of strat. water vapour fix                           MORC.387   
           ENDIF                                                             MORC.388   
                                                                             MORC.389   
c Ozone                                                                      MORC.390   
                                                                             MORC.391   
           o3(l)=AMFRAC*o3mod2(ld,i,ihem,JH)+                                MORC.392   
     +            (1.0-AMFRAC)*o3mod1(ld,i,ihem,JH)                          MORC.393   
                                                                             MORC.394   
                                                                             MORC.395   
          ENDDO                             ! End of loop over column.       MORC.396   
c ----------------------------------------------------- And alat1            MORC.397   
                                                                             MORC.398   
          alat1=alat(JH)*REAL(-(ihem*2.)+3)                                  MORC.399   
                                                                             MORC.400   
c 9-7-97 Piers' clouds                                                       MORC.401   
c                                                                            MORC.402   
          cf(1,1)=max(cfrac(im,5),cfrac(im,4))  !deep and shallow            MORC.403   
          cf(2,1)=cfrac(im,3)    !high cloud fraction                        MORC.404   
          cf(3,1)=cfrac(im,2)    ! mid cloud fraction                        MORC.405   
          cf(4,1)=cfrac(im,1)    !low cloud fraction                         MORC.406   
                                                                             MORC.407   
C cloud levels reversed                                                      MORC.408   
c choose lowest level for bottom                                             MORC.409   
          ic(1,1)=max(icflag(im,4,1),icflag(im,5,1)) !bottom deep/shallow    MORC.410   
C choose hightest level for top                                              MORC.411   
          ic(1,2)=min(icflag(im,4,2),icflag(im,5,2)) !top deep/shallow       MORC.412   
C reverse                                                                    MORC.413   
          ic(1,1)=NLP-ic(1,1)+1                                              MORC.414   
          ic(1,2)=NLP-ic(1,2)+1                                              MORC.415   
C makes sure bottom cloud isn't in bottom level                              MORC.416   
          if (ic(1,1).eq.2) ic(1,1)=3                                        MORC.417   
          if (ic(1,2).eq.2) ic(1,1)=3                                        MORC.418   
                                                                             MORC.419   
          ic(2,1)=NLP-icflag(im,3,1)+1 ! high pos, one height needed         MORC.420   
          ic(3,1)=NLP-icflag(im,2,1)+1 ! mid pos                             MORC.421   
          ic(4,1)=NLP-icflag(im,1,1)+1 ! low pos                             MORC.422   
                                                                             MORC.423   
          SWALB=SALB(IM,JH)                                                  MORC.424   
C cloud cf and ic passed. fluxes returned.                                   MORC.425   
C which is net flux at TOA in profile                                        MORC.426   
C Call radiation scheme                                                      MORC.427   
          alon=REAL(i-1)/REAL(mg)*360.0                                      MORC.428   
          call nikosrad(pr,t,h2o,o3,alat1,htlw,htsw,DOY,cf,ic,               MORC.429   
     $     fluxes,swalb,alon)                                                MORC.430   
                                                                             MORC.431   
c store net flux in PNET                                                     MORC.432   
          PNET(IM,JH)=fluxes(1,1,1)-fluxes(1,2,1)+fluxes(2,1,1)-             MORC.433   
     $     fluxes(2,2,1)                                                     MORC.434   
          SNET(IM,JH)=fluxes(1,1,2)-fluxes(1,2,2)+fluxes(2,1,2)-             MORC.435   
     $     fluxes(2,2,2)                                                     MORC.436   
          rrflux(im,jh,1)=fluxes(1,1,2)                                      MORC.437   
          rrflux(im,jh,2)=fluxes(1,2,2)                                      MORC.438   
          rrflux(im,jh,3)=fluxes(2,1,2)                                      MORC.439   
          rrflux(im,jh,4)=fluxes(2,2,2)                                      MORC.440   
          rrflux(im,jh,5)=fluxes(1,1,1)-fluxes(1,2,1)                        MORC.441   
          rrflux(im,jh,6)=fluxes(2,2,1)                                      MORC.442   
                                                                             MORC.443   
          DO l=nl,1,-1                                                       MORC.444   
c  bottom heating rate is zero in morecret                                   MORC.445   
           LD=NL+1-L                                                         MORC.446   
           IM=I+IOFM                                                         MORC.447   
           HTNETO=HTNET(IHem,JH,I,LD)                                        MORC.448   
           htnet(ihem,jh,i,ld)=htlw(l+1)+htsw(l+1)                           MORC.449   
                                                                             MORC.450   
c sets this heating rate                                                     MORC.451   
           TTRD(IM,LD)=(HTNETO                                               MORC.452   
     $      +HTNET(IHEM,JH,I,LD))/(CHRF*2.)                                  MORC.453   
                                                                             MORC.454   
c  put in linear interpolation of heating rates between this                 MORC.455   
c  longitude and last one calculated (i-nskip)                               MORC.456   
           IF (i-ilast.gt.1) then                                            MORC.457   
            DO j=ilast+1,i-1                                                 MORC.458   
              a=REAL(j-ilast)/REAL(i-ilast)                                  MORC.459   
              b=1.-a                                                         MORC.460   
              HTNETO=HTNET(IHEM,JH,J,LD)                                     MORC.461   
              htnet(ihem,jh,j,ld)=a*htnet(ihem,jh,i,ld)+                     MORC.462   
     $                         b*htnet(ihem,jh,ilast,ld)                     MORC.463   
              im=j+iofm                                                      MORC.464   
              TTRD(IM,LD)=(HTNETO                                            MORC.465   
     $                   +HTNET(IHEM,JH,J,LD))/(CHRF*2.)                     MORC.466   
              IF (l.eq.nl) then                                              MORC.467   
                pnet(im,jh)=a*pnet(i+iofm,jh)+                               MORC.468   
     $                     b*pnet(ilast+iofm,jh)                             MORC.469   
                snet(im,jh)=a*snet(i+iofm,jh)+                               MORC.470   
     $                     b*snet(ilast+iofm,jh)                             MORC.471   
                DO k=1,6                                                     MORC.472   
                  rrflux(im,jh,k)=a*rrflux(i+iofm,jh,k)                      MORC.473   
     $                           +b*rrflux(ilast+iofm,jh,k)                  MORC.474   
                ENDDO                                                        MORC.475   
              ENDIF                                                          MORC.476   
            ENDDO                                                            MORC.477   
           ENDIF                                                             MORC.478   
          ENDDO                                                              MORC.479   
          ilast=i                                                            MORC.480   
C end of conditional execution of morcrette code                             MORC.481   
         ENDIF                                                               MORC.482   
C end of loop over longitudes                                                MORC.483   
        ENDDO                                                                MORC.484   
        IF (ilast.ne.mg) then                                                MORC.485   
         DO j=ilast+1,mg                                                     MORC.486   
           a=REAL(j-ilast)/REAL(mg+1-ilast)                                  MORC.487   
           b=1.-a                                                            MORC.488   
           im=j+iofm                                                         MORC.489   
           DO l=nl,1,-1                                                      MORC.490   
             ld=nl+1-l                                                       MORC.491   
             HTNETO=HTNET(IHEM,JH,J,LD)                                      MORC.492   
             htnet(ihem,jh,j,ld)=a*htnet(ihem,jh,1,ld)+                      MORC.493   
     $                         b*htnet(ihem,jh,ilast,ld)                     MORC.494   
             TTRD(IM,LD)=(HTNET(IHEM,JH,J,LD)                                MORC.495   
     $                    +HTNETO)/(CHRF*2.)                                 MORC.496   
             IF (l.eq.nl) then                                               MORC.497   
               pnet(im,jh)=a*pnet(1+iofm,jh)+                                MORC.498   
     $                     b*pnet(ilast+iofm,jh)                             MORC.499   
               snet(im,jh)=a*snet(1+iofm,jh)+                                MORC.500   
     $                     b*snet(ilast+iofm,jh)                             MORC.501   
               DO k=1,6                                                      MORC.502   
                rrflux(im,jh,k)=a*rrflux(1+iofm,jh,k)                        MORC.503   
     $                          +b*rrflux(ilast+iofm,jh,k)                   MORC.504   
               ENDDO                                                         MORC.505   
             ENDIF                                                           MORC.506   
           ENDDO                                                             MORC.507   
         ENDDO                                                               MORC.508   
        ENDIF                                                                MORC.509   
c       if (ihem.eq.2) print *, 'rad ',jh,(pnet(im1,jh),im1=1,IGC)           MORC.510   
       ELSE                                                                  MORC.511   
C                                                                            MORC.512   
c  Doesn't do rad scheme (simply uses old heating rates)                     MORC.513   
C                                                                            MORC.514   
         DO i=1,mg                                                           MORC.515   
           DO LD=1,NL                                                        MORC.516   
             im=i+IOFM                                                       MORC.517   
             TTRD(im,LD)=(htnet(ihem,jh,i,ld))/CHRF                          MORC.518   
           ENDDO                                                             MORC.519   
         ENDDO                                                               MORC.520   
       ENDIF                                                                 MORC.521   
       IOFM=MGPP                                                             MORC.522   
 800  CONTINUE ! end of loop over hemispheres                                MORC.523   
      IF (LSHORT.AND.(KOUNT.eq.1)) then                                      MORC.524   
        DO l=1,nl                                                            MORC.525   
          DO i=1,igc                                                         MORC.526   
            ttrd(i,l)=ttrd(i,l)*2.                                           MORC.527   
          ENDDO                                                              MORC.528   
        ENDDO                                                                MORC.529   
      ENDIF                                                                  MORC.530   
                                                                             MORC.531   
      RETURN                                                                 MORC.532   
      END                                                                    MORC.533   
