program mkpham_gswm09 C pgf90 -o mkpham_gswm09.exe mkpham_gswm09.f C 05/13 bae: Make monthly gswm09 files of D(iurnal) and S(emidiurnal) C values amp+ph of geopotential ht, Tn, Un and Vn on sres/dres glat/glon 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 1_month of amp(ph)(lat,lon,ht,4or6) 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 :: ampk(144,72,13,6),phk(144,72,13,6),fac1(72),fak1(13) real :: glat09(59),ht09(55),glon09(72),temp(6) real :: glattie(72),glontie(144),ht2p5(13),glon144(144) real :: ampjk1,ampjk2,phjk1,phjk2 real :: adphk0,adphk1,adphk10,adphk11,adphk20,adphk21 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*30 :: outfile character*22 :: file95km character*1 :: blank integer :: jj,jr real :: utx 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' do k=1,13 ht2p5(k) = 90. + (k-1)*2.5 enddo 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 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_dresD_90-120km_pham_',a3)") mon write (file95km,"('gswm09_D_95km_pham_',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_dresS_90-120km_pham_',a3)") mon write (file95km,"('gswm09_S_95km_pham_',a3)") mon endif open (9,file=file95km) 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 C 9/14: Not necessary to Skip over blank after each ht because of * for 6 values C read (ird,"(a1)") blank C read (ird+1,"(a1)") blank C write (6,"(1x,'ird k temp(6)=',2i3,e12.4)") ird,k,temp(6) enddo enddo ! do ird=10,20,2 C 9/14: Print out amps and phases for GSMW09 grid to compare with TIE grid ratty amp/ph k95 = 24 if (nn==1) then write (9,"(f8.2,' km GSWM09 D for 72 glons and 59', | ' glats in ',a3)") ht09(k95),mon else write (9,"(f8.2,' km GSWM09 SD for 72 glons and 59', | ' glats in ',a3)") ht09(k95),mon endif write (9,"(10f8.2)") (glon09(i),i=37,72),(glon09(j),j=1,36) write (9,"(10f8.2)") glat09 C Print out 95km ph and amp of zg in km, Tn in K, and Un,Vn in m/s write (9,"(' ph amp: gph(km) tn(K) un_vn (m/s)')") do i=37,72 do j=1,59 write (9,"(4(f6.2,e13.5))") | ph(i,j,k95,1),amp(i,j,k95,1)*0.001, | ph(i,j,k95,3),amp(i,j,k95,3), | ph(i,j,k95,4),amp(i,j,k95,4), | ph(i,j,k95,5),amp(i,j,k95,5) enddo enddo do i=1,36 do j=1,59 write (9,"(4(f6.2,e13.5))") | ph(i,j,k95,1),amp(i,j,k95,1)*0.001, | ph(i,j,k95,3),amp(i,j,k95,3), | ph(i,j,k95,4),amp(i,j,k95,4), | ph(i,j,k95,5),amp(i,j,k95,5) enddo enddo C 9/14: revise ph=0 when amp=0 near poles to be same ph as close glats so can interpolate correctly do i=1,72 do k=1,55 do n=1,6 C Do not extrapolate, just set to same as first non-zero value C Near SH do jr=1,4 j = 5-jr if (ph(i,j,k,n) .lt. 0.00001) then do jj=1,j ph(i,jj,k,n) = ph(i,j+1,k,n) enddo endif enddo C Near NH do j=56,59 if (ph(i,j,k,n) .lt. 0.00001) then do jj=j,59 ph(i,jj,k,n) = ph(i,j-1,k,n) enddo endif enddo enddo enddo enddo C 9/14: Print out amps and phases for GSMW09 grid to compare with TIE grid ratty amp/ph k95 = 24 if (nn==1) then write (9,"(f8.2,' km GSWM09 D for 72 glons and 59', | ' glats in ',a3)") ht09(k95),mon else write (9,"(f8.2,' km GSWM09 SD for 72 glons and 59', | ' glats in ',a3)") ht09(k95),mon endif write (9,"(10f8.2)") (glon09(i),i=37,72),(glon09(j),j=1,36) write (9,"(10f8.2)") glat09 C Print out 95km ph and amp of zg in km, Tn in K, and Un,Vn in m/s write (9,"(' ph amp: gph(km) tn(K) un_vn (m/s)')") do i=37,72 do j=1,59 write (9,"(4(f6.2,e13.5))") | ph(i,j,k95,1),amp(i,j,k95,1)*0.001, | ph(i,j,k95,3),amp(i,j,k95,3), | ph(i,j,k95,4),amp(i,j,k95,4), | ph(i,j,k95,5),amp(i,j,k95,5) enddo enddo do i=1,36 do j=1,59 write (9,"(4(f6.2,e13.5))") | ph(i,j,k95,1),amp(i,j,k95,1)*0.001, | ph(i,j,k95,3),amp(i,j,k95,3), | ph(i,j,k95,4),amp(i,j,k95,4), | ph(i,j,k95,5),amp(i,j,k95,5) enddo enddo C real :: ampk(144,72,13,6),phk(144,72,13,6) C Print out the interpolated amps and phases for zg,tn,un,vn at dres grid utx = 24. if (nn .eq. 2) utx = 12. 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 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 C ampjk1 = fac1(j)*amp(i2,j1(j),k1(k),n) + C | (1.-fac1(j))*amp(i2,j1(j),k1(k)+1,n) C ampjk2 = fac1(j)*amp(i2,jp1,k1(k),n) + C | (1.-fac1(j))*amp(i2,jp1,k1(k)+1,n) ampjk1 = max(0.,fac1(j)*amp(i2,j1(j),k1(k),n) + | (1.-fac1(j))*amp(i2,jp1,k1(k),n)) ampjk2 = max(0.,fac1(j)*amp(i2,j1(j),k1(k)+1,n) + | (1.-fac1(j))*amp(i2,jp1,k1(k)+1,n)) ampk(i,j,k,n) = fak1(k)*ampjk1 + (1.-fak1(k))*ampjk2 ! Check to see if phase crosses 0/24h for nn=1, or 0/12h for nn=2 (added 9/14) adphk10=0. C if (ph(i2,j1(j),k1(k),n)-ph(i2,j1(j),k1(k)+1,n)<-12.)adphk10=24 if (ph(i2,j1(j),k1(k),n)-ph(i2,jp1,k1(k),n)<-utx/2.)adphk10=utx adphk11=0. C if (ph(i2,j1(j),k1(k),n)-ph(i2,j1(j),k1(k)+1,n)>12.)adphk11=24 if (ph(i2,j1(j),k1(k),n)-ph(i2,jp1,k1(k),n)>utx/2.)adphk11=utx C phjk1 = fac1(j)*(ph(i2,j1(j),k1(k),n)+adphk10) + C | (1.-fac1(j))*(ph(i2,j1(j),k1(k)+1,n)+adphk11) phjk1 = fac1(j)*(ph(i2,j1(j),k1(k),n)+adphk10) + | (1.-fac1(j))*(ph(i2,jp1,k1(k),n)+adphk11) adphk20=0. C if (ph(i2,jp1,k1(k),n)-ph(i2,jp1,k1(k)+1,n)<-12.)adphk20=24 if (ph(i2,j1(j),k1(k)+1,n)-ph(i2,jp1,k1(k)+1,n)<-utx/2.) | adphk20=utx adphk21=0. C if (ph(i2,jp1,k1(k),n)-ph(i2,jp1,k1(k)+1,n)>12.)adphk21=24 if (ph(i2,j1(j),k1(k)+1,n)-ph(i2,jp1,k1(k)+1,n)>utx/2.) | adphk21=utx C phjk2 = fac1(j)*(ph(i2,jp1,k1(k),n)+adphk20) + C | (1.-fac1(j))*(ph(i2,jp1,k1(k)+1,n)+adphk21) phjk2 = fac1(j)*(ph(i2,j1(j),k1(k)+1,n)+adphk20) + | (1.-fac1(j))*(ph(i2,jp1,k1(k)+1,n)+adphk21) adphk0=0. if (phjk1-phjk2<-utx/2.) adphk0 = utx adphk1=0. if (phjk1-phjk2>utx/2.) adphk1 = utx phk(i,j,k,n) = fak1(k)*(phjk1+adphk0) + (1.-fak1(k))* | (phjk2+adphk1) ! Make phase 0-utxh if (phk(i,j,k,n)<0.) phk(i,j,k,n)=phk(i,j,k,n)+utx if (phk(i,j,k,n)>=utx) phk(i,j,k,n)=phk(i,j,k,n)-utx 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 ampk(i,j,k,n) = 0.5*(ampk(i-1,j,k,n) + ampk(ip1,j,k,n)) adphk0=0. if (phk(i-1,j,k,n)-phk(ip1,j,k,n)<-utx/2.) adphk0 = utx adphk1=0. if (phk(i-1,j,k,n)-phk(ip1,j,k,n)>utx/2.) adphk1 = utx phk(i,j,k,n) = 0.5*(phk(i-1,j,k,n)+adphk0 + phk(ip1,j,k,n) | +adphk1) ! Make phase 0-utx h if (phk(i,j,k,n)<0.) phk(i,j,k,n)=phk(i,j,k,n)+utx if (phk(i,j,k,n)>=utx) phk(i,j,k,n)=phk(i,j,k,n)-utx enddo enddo enddo enddo C Print out outfile ph and amp zg in cm (*100), Tn in K, and Un,Vn in cm/s (*100) do i=1,144 do k=1,13 do j=1,72 write (30,"(4(f6.2,e13.5))") | phk(i,j,k,1),ampk(i,j,k,1)*100., | phk(i,j,k,3),ampk(i,j,k,3), | phk(i,j,k,4),ampk(i,j,k,4)*100., | phk(i,j,k,5),ampk(i,j,k,5)*100. enddo enddo enddo enddo ! do nn=1,2 stop end