program mkgswm09

C pgf90 -o mkgswm09.exe mkgswm09.f
C 04/13 bae:  Make monthly gswm09 files of D(iurnal) and S(emidiurnal)
C  values vs UT/glon of geopotential ht, relative density, Tn,
C  and neutral winds Un, Vn and Wn on the tiegcm double resolution grid
C  from 90 to 120 km at 13 hts 2.5km apart.  Do not use the glat/glon
C  grid of Zhang et al. [2010a,b] JGR or their ht grid.

C  Phase for D ~12LT (peak at noon), and for S ~0 or 12LT (same)
C   Save as (lon,lat,1month,UThr)
C   value = amp*cos(nn*pi*(LT-phase)/12), where nn=1D,2S for each UT/glon

      real :: amp(72,59,55,6),ph(72,59,55,6),xamp(6)
      real :: valtie(144,72,13,6,24),fac1(72),fak1(13)
      real :: glat09(59),ht09(55),glon09(72),temp(6)
      real :: glattie(72),glontie(144),ht2p5(13),slt1,pi,glon144(144)
      real :: valj1k1,valj1k2,valj2k1,valj2k2,valjk1,valjk2
      integer :: i,j,k,m,n,nn,iut,j1(72),k1(13),jp1,jx(6),ix(6)
      character*3 :: mon,mon12(12)
      character*22 :: flzgda,flrrda,fltnda,flunda,flvnda,flwnda
      character*22 :: flzgdp,flrrdp,fltndp,flundp,flvndp,flwndp
      character*22 :: flzgsa,flrrsa,fltnsa,flunsa,flvnsa,flwnsa
      character*22 :: flzgsp,flrrsp,fltnsp,flunsp,flvnsp,flwnsp
      character*29 :: outfile
      character*42 :: file95km
      character*12 :: monfile
      data mon12/'jan','feb','mar','apr','may','jun','jul','aug',
     |           'sep','oct','nov','dec'/
C  #22-30=86.293 to 120.019km so interpolate to 90:2.5:120 from these
C #24=94.58km
      data ht09/0.,4.161,8.188,12.201,16.218,20.239,24.264,28.295,
     |  32.332,36.378,40.434,44.503,48.588,52.691,56.819,60.975,65.167,
     |  69.403,73.695,78.056,82.129,86.293,90.57,94.58,98.736,
     |  103.075,107.177,111.514,115.619,120.019,124.175,128.688,
     |  132.911,137.539,141.772,146.414,150.470,154.871,159.653,
     |  164.847,169.021,173.437,178.086,182.957,188.030,193.283,
     |  198.689,204.224,209.865,215.589,221.379,227.221,233.102,
     |  239.013,244.946/
      
     
      pi = 2.*asin(1.)
      do k=1,13
       ht2p5(k) = 90. + (k-1)*2.5
      enddo
C      do i=1,44	! error found 9/14 and glon09 180 off from glontie!!
C	glontie(i) = -180. + (i-1)*2.5
      do i=1,144
C 9/14:  12h off because glontie is 180 deg off from glon09 so i2=(i+1)/2 is wrong!
	glontie(i) = -180. + (i-1)*2.5
        glon144(i) = glontie(i)
        if (glon144(i) .lt. 0.) glon144(i) = glon144(i) + 360.
      enddo
      do i=1,72
        glon09(i) = (i-1)*5
        glattie(i) = -88.75 + (i-1)*2.5
      enddo
      do j=1,59
        glat09(j) = -90 + j*3
      enddo

C Find j1 and fac1 for interpolation in glat from glat09 to glattie
C  Must extrapolate at j=1,72
      j1(1) = 1
      fac1(1) = (abs(glattie(1))-abs(glat09(2)))/3.
      j1(72) = 59
      fac1(72) = fac1(1) 
      do j=2,36
        j2 = 72 - j + 1
        do k=1,29
	 if (abs(glattie(j)).le.abs(glat09(k)) .and. abs(glattie(j))
     |  .ge. abs(glat09(k+1))) go to 100
        enddo
        write (6,"(1x,'stop j,glattie(j),abs(glat09(1,30))=',i2,3f7.2 
     |   )") j,glattie(j),glat09(1),glat09(30)
        stop
100     j1(j) = k
        j1(j2) = 59-k+1
        fac1(j) = (abs(glattie(j))-abs(glat09(k+1)))/3.
        fac1(j2) = fac1(j)
      enddo

C Find k1 and fak1 for interpolation in ht from ht09 to ht2p5
      do k=1,13
        do j=22,29
          if(ht2p5(k).ge.ht09(j) .and. ht2p5(k).le.ht09(j+1)) go to 101
        enddo
        write (6,"(1x,'stop k,ht2p5(k),ht09(22),ht09(30)=',i2,3f7.2)") 
     |   k,ht2p5(k),ht09(22),ht09(30)
        stop
101     k1(k) = j
        fak1(k) = (ht09(j+1)-ht2p5(k))/(ht09(j+1)-ht09(j))
      enddo

      k95 = 3

C TEMP: Do only 1 month
      do mond=3,3
      mon = mon12(mond)
      write (monfile,"('max_amps_',a3)") mon
      open (31,file=monfile)

      write (file95km,"('gswm09_gph_relrho_tn_un_vn_wn_95km_0ut_',a3)")
     |  mon
      open (9,file=file95km)
      write (9,"(f8.2,' km GSWM09 D then SD for 144 glons and 72', 
     |  ' glats in ',a3)") ht2p5(k95),mon
      write (9,"(10f8.2)") glontie
      write (9,"(10f8.2)") glattie
C Print out 95km zg in m, rr in %, Tn in K, and Un,Vn,Wn in m/s 
      write (9,"(' pert gph(m) relrho(%) tn(K) un_vn_wn (m/s)')")

C Open 24 kinds of files for 'mon'
      write (flzgda,"('geoptt_',a3,'_all_Dam.dat')") mon
      write (flzgdp,"('geoptt_',a3,'_all_Dph.dat')") mon
      write (flzgsa,"('geoptt_',a3,'_all_Sam.dat')") mon
      write (flzgsp,"('geoptt_',a3,'_all_Sph.dat')") mon
      write (flrrda,"('rlrrho_',a3,'_all_Dam.dat')") mon
      write (flrrdp,"('rlrrho_',a3,'_all_Dph.dat')") mon
      write (flrrsa,"('rlrrho_',a3,'_all_Sam.dat')") mon
      write (flrrsp,"('rlrrho_',a3,'_all_Sph.dat')") mon
      write (fltnda,"('temper_',a3,'_all_Dam.dat')") mon
      write (fltndp,"('temper_',a3,'_all_Dph.dat')") mon
      write (fltnsa,"('temper_',a3,'_all_Sam.dat')") mon
      write (fltnsp,"('temper_',a3,'_all_Sph.dat')") mon
      write (flunda,"('zonwin_',a3,'_all_Dam.dat')") mon
      write (flundp,"('zonwin_',a3,'_all_Dph.dat')") mon
      write (flunsa,"('zonwin_',a3,'_all_Sam.dat')") mon
      write (flunsp,"('zonwin_',a3,'_all_Sph.dat')") mon
      write (flvnda,"('merwin_',a3,'_all_Dam.dat')") mon
      write (flvndp,"('merwin_',a3,'_all_Dph.dat')") mon
      write (flvnsa,"('merwin_',a3,'_all_Sam.dat')") mon
      write (flvnsp,"('merwin_',a3,'_all_Sph.dat')") mon
      write (flwnda,"('verwin_',a3,'_all_Dam.dat')") mon
      write (flwndp,"('verwin_',a3,'_all_Dph.dat')") mon
      write (flwnsa,"('verwin_',a3,'_all_Sam.dat')") mon
      write (flwnsp,"('verwin_',a3,'_all_Sph.dat')") mon

C  Make 2 files:  diurnal and semidiurnal
      do nn=1,2
       write (31,"(1x,'nn =',i1)") nn
       if (nn==1) then
        open (10,file=flzgda,status='old')
        open (11,file=flzgdp,status='old')
        open (12,file=flrrda,status='old')
        open (13,file=flrrdp,status='old')
        open (14,file=fltnda,status='old')
        open (15,file=fltndp,status='old')
        open (16,file=flunda,status='old')
        open (17,file=flundp,status='old')
        open (18,file=flvnda,status='old')
        open (19,file=flvndp,status='old')
        open (20,file=flwnda,status='old')
        open (21,file=flwndp,status='old')
        write (outfile,"('gswm09_dres_90-120km_diur_',a3)") mon
       else
        open (10,file=flzgsa,status='old')
        open (11,file=flzgsp,status='old')
        open (12,file=flrrsa,status='old')
        open (13,file=flrrsp,status='old')
        open (14,file=fltnsa,status='old')
        open (15,file=fltnsp,status='old')
        open (16,file=flunsa,status='old')
        open (17,file=flunsp,status='old')
        open (18,file=flvnsa,status='old')
        open (19,file=flvnsp,status='old')
        open (20,file=flwnsa,status='old')
        open (21,file=flwnsp,status='old')
        write (outfile,"('gswm09_dres_90-120km_semi_',a3)") mon
       endif
       open (30,file=outfile)

      n = 0
      do ird=10,20,2
       n = n + 1
      do k=1,55
      do j=1,59
       do i=1,72,6
         read (ird,*) temp
          do m=1,6
           amp(i+m-1,j,k,n) = temp(m)
          enddo
         read (ird+1,*) temp
          do m=1,6
           ph(i+m-1,j,k,n) = temp(m)
          enddo
       enddo
      enddo
      enddo
      enddo	! do ird=10,20,2

C Print out max amp for each in hts 22-30
      do k=22,30
       do n=1,6
        xamp(n) = 0.
         do i=1,72
          do j=1,59
           if (xamp(n).lt.amp(i,j,k,n)) then
            xamp(n)=amp(i,j,k,n)
            jx(n) = j
            ix(n) = i
           endif
          enddo
         enddo
       enddo
       write (31,"(1x,'km xamp are:',f7.2,6e13.5/
     |   '  at jx ix =',6i3,1x,6i4)") ht09(k),xamp,jx,ix
      enddo


C   real :: valtie(144,72,13,24,6)
C   value = amp*cos(nn*pi*(LT-phase)/12), where nn=1D,2S for each UT/glon
C  Calculate the value for 24 UTs at each glon position

      do iut=0,23
       do i=1,143,2
C 9/14:  12h off because glontie is 180 deg off from glon09 so i2=(i+1)/2 is wrong!
         i2 = (i+1)/2
	 i2 = i2 + 36
         if (i2 .gt. 72) i2 = i2 - 72
         if (abs(glon09(i2)-glon144(i)) .gt. 0.01) then
          write (6,"(1x,'bad glon match for glon09(i2) glon144(i) =',
     |     2f8.2)") glon09(i2),glon144(i)
	  stop
	 endif
         slt1 = iut + glontie(i)/15.
         if (slt1.lt.0.) slt1 = slt1 + 24.
C         write (6,"(1x,'iut glon slt1 =',i2,2f8.2)")iut,glontie(i),slt1
C TEMP 9/14: Substitute UT for LT since UT is used in gswm09.F for TIEGCM (???)
C  LT makes D look like SD, and SD like 4D, while UT looks better
         slt1 = iut
         do j=1,72
C SH
          jp1 = j1(j) + 1
C NH
          if (j .gt. 36) jp1 = j1(j) - 1
          do k=1,13
           do n=1,6
            valj1k1 = amp(i2,j1(j),k1(k),n)*cos(nn*pi*(slt1-
     |       ph(i2,j1(j),k1(k),n))/12.)
            valj1k2 = amp(i2,j1(j),k1(k)+1,n)*cos(nn*pi*(slt1-
     |       ph(i2,j1(j),k1(k)+1,n))/12.)
            valj2k1 = amp(i2,jp1,k1(k),n)*cos(nn*pi*(slt1-
     |       ph(i2,jp1,k1(k),n))/12.)
            valj2k2 = amp(i2,jp1,k1(k)+1,n)*cos(nn*pi*(slt1-
     |       ph(i2,jp1,k1(k)+1,n))/12.)
            valjk1 = fac1(j)*valj1k1 + (1.-fac1(j))*valj2k1
            valjk2 = fac1(j)*valj1k12+ (1.-fac1(j))*valj2k2
            valtie(i,j,k,n,iut+1) = fak1(k)*valjk1+(1.-fak1(k))*valjk2
           enddo
          enddo
         enddo
       enddo
C Interpolate between glons
       do i=2,144,2
         ip1 = i + 1
         if (i==144) ip1 = 1
         do j=1,72
          do k=1,13
           do n=1,6
            valtie(i,j,k,n,iut+1) = 0.5*(valtie(i-1,j,k,n,iut+1) +
     |       valtie(ip1,j,k,n,iut+1))
           enddo
          enddo
         enddo
       enddo
C Print out outfile zg in cm (*100), rr in %, Tn in K, and Un,Vn,Wn in cm/s (*100) 
       do i=1,144
         do j=1,72
          if (iut==0) 
     |      write (9,"(6e13.5)") valtie(i,j,k95,1,iut+1),
     |      valtie(i,j,k95,2,iut+1),valtie(i,j,k95,3,iut+1),
     |      valtie(i,j,k95,4,iut+1),valtie(i,j,k95,5,iut+1),
     |      valtie(i,j,k95,6,iut+1)
          do k=1,13
           write (30,"(6e13.5)") valtie(i,j,k,1,iut+1)*100.,
     |      valtie(i,j,k,2,iut+1),valtie(i,j,k,3,iut+1),
     |      valtie(i,j,k,4,iut+1)*100.,valtie(i,j,k,5,iut+1)*100.,
     |      valtie(i,j,k,6,iut+1)*100.
          enddo
         enddo
       enddo
       write (6,"(1x,'end iut =',i2)") iut
      enddo	! do iut=0,23

C  Find max vals and jx,ix
      do k=1,13
       do n=1,6
        xamp(n) = 0.
       do iut=1,24
       do i=1,143,2
        do j=1,72
           if (xamp(n).lt.valtie(i,j,k,n,iut)) then
            xamp(n)=valtie(i,j,k,n,iut)
            jx(n) = j
            ix(n) = i
           endif
        enddo
       enddo
       enddo
       enddo
       write (31,"(1x,'tie km xval are:',f7.2,6e13.5/
     |   '  at jx ix =',6i3,1x,6i4)") ht2p5(k),xamp,jx,ix
      enddo        
    
      close (30)
      enddo	! do nn=1,2
      close (31)
      enddo	! do mond=1,12

      stop
      end