program mkgswm09 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 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 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 C #22-30=86.293 to 120.019km so interpolate to 90:2.5:120 from these 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/ mon = 'mar' open (31,file='max_amps_mar') pi = 2.*asin(1.) do k=1,13 ht2p5(k) = 90. + (k-1)*2.5 enddo do i=1,144 glontie(i) = -180. + (i-1)*2.5 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 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 i2 = (i+1)/2 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 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 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 enddo ! do nn=1,2 stop end