前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >气象编程 | 大气视热源Q1和视水汽汇Q2程序

气象编程 | 大气视热源Q1和视水汽汇Q2程序

作者头像
bugsuse
发布2020-12-16 11:52:01
2.9K0
发布2020-12-16 11:52:01
举报
文章被收录于专栏:气象杂货铺气象杂货铺
大气视热源Q1和视水汽汇Q2的程序

提供NCL、FORTRAN、IDL三种选择,推荐使用NCL。

Q1Q2_yanai.ncl: Using high-frequency [eg: 3/6/12 hourly or daily] data, calculate apparent-heat-source (Q1) and apparent-moisture-sink (Q2) quantities as described in:

This example uses daily mean data from NOAA/OAR/ESRL PSD, Boulder, Colorado, USA. Specifically, NCEP_Reanalysis 2 daily mean data spanning 3-9 April 1993.

This is a test script. It is not fully tested. The returned quantities are:

代码语言:javascript
复制
         q1 =   (dT/dt) - [omega*static_stability - V.grad(T)]  ; K/day
         q2 = -[(dH/dt) +  V.grad(H) +omega*(dH/dp)]            ; g/(kg-Pa)
    and
         Q1 = Cp*q1     ; apparent Heat Source
         Q2 = Lc*q2     ; apparent Moisture Sink     

具体NCL代码如下:

代码语言:javascript
复制

undef("Q1Q2_yanai.ncl")
function Q1Q2_yanai(time[*]:numeric,p,u,v,H,T,omega,npr[1]:integer,opt[1]:logical)
;+++++++++++++++++++++++++++++++++++++++++++++++++++++
; NOT SUPPORTED:  NOT FULLY TESTED : WORK in PROGRESS
;                  **CHECK UNITS**
;+++++++++++++++++++++++++++++++++++++++++++++++++++++
; Nomenclature
; time    - "seconds since ..."
; p       - Pressure [Pa]
; u,v     - zonal, meridional wind components[m/s]
; H       - specific humidity [g/kg]
; T       - temperature [K]  or [C]
; omega   - vertical velocity [Pa/s]
; npr     - demension number corresponding to pressure dimension
; opt     - set to False
;+++++++++++++++++++++++++++++++++++++++++++++++++++++
;;  Q1  = Cp*dTdt - Cp*(omega*ss-advT)
;;  Q1  = Cp*[ dTdt- omega*ss-advT ]
;;  q1  = dTdt- omega*ss-advT    
;;  Q1  = Cp*[ q1 ]                   
;;
;;  Q1  = Total diabatic heating [including radiation, latent heating, and
;;        sfc heat flux) and sub-grid scale heat flux convergence
;;---
;;  q2  = -(dHdt +advH +dHdp)
;;  Q2  = Lc*[ q2 ]
;;
;;  Q2  = the latent heating due to condensation or evaporation processes
;;        and subgrid-scale moisture flux convergences,
;++++++++++++++++++++++++++++++++++++++++++++
; References:
; 
; https://renqlsysu.github.io/2019/02/01/apparent_heat_source/
;
; Yanai, M. (1961): 
; A detailed analysis of typhoon formation.
; J. Meteor. Soc. Japan , 39 , 187–214
; 
; Yanai et al (1973):
;  Determination of bulk properties of tropical cloud clusters 
;    from large-scale heat and moisture budgets.
;  J.  Atmos.  Sci.  , 30 , 611–627,
;  https://doi.org/10.1175/1520-0469(1973)030<0611:DOBPOT>2.0.CO;2
;
; Yanai, M and T.Tomita (1998):
; https://pdfs.semanticscholar.org/fb57/a6a59cc4a684194b5e622ea6f875d0b4439a.pdf
;********************************************
; Diabatic heating in the atmosphere is a combined consequence of 
; radiative fluxes, phase changes of water substance, and turbulent 
; flux of sensible heat from the earth's surface. 
; 
; In the tropics, it is the major driving force of the atmospheric circulation. 
; It responds to the vertical gradient of diabatic heating.
;********************************************
local dimp, dimu, dimv, dimH, dimT, dimo             \
    , rankp, ranku, rankv, rankH, rankT, ranko       \
    , Cp, Lc, dTdt, ss, opt_adv, long_name, units    \
    , gridType, advT, q1, Q1, dHdt, advH, loq, q2, Q2
begin
;;             Use for error testing
;;dimp         = dimsizes(p)
;;dimu         = dimsizes(u)
;;dimv         = dimsizes(v)
;;dimH         = dimsizes(H)
;;dimT         = dimsizes(T)
;;dimo         = dimsizes(omega)

;;rankp        = dimsizes(dimp)
;;ranku        = dimsizes(dimu)
;;rankv        = dimsizes(dimv)
;;rankH        = dimsizes(dimH)
;;rankT        = dimsizes(dimT)
;;ranko        = dimsizes(dimo)

;*******************************************
;---Compute local dT/dt  [K/s]
;*******************************************

  dTdt = center_finite_diff_n (T,time,False,0,0)   ; 'time' is 'seconds since'
  copy_VarCoords(T, dTdt)
  dTdt@longname = "Temperature: Local Time derivative"
  dTdt@units = "K/s"
  printVarSummary(dTdt)
  printMinMax(dTdt,0)
  print("-----")
   
;****************************************
;---Compute static stability [K/Pa] <==>  or "K-m-s2/kg"
;****************************************

  ss  = static_stability (p , T, 1, 0) 
  printVarSummary(ss)
  printMinMax(ss,0)
  print("-----")

;****************************************
;---Compute advection term: spherical harmonics
;---U*(dT/dx) + V*(dT/dy):  [m/s][K/m] -> [K/s]
;****************************************

  opt_adv   = 0
  long_name = "temp advection"
  units     = "K/s"
  gridType  = 1  
  advT      = advect_variable(u,v,T,gridType,long_name,units,opt_adv) 
   
;****************************************
;---Net Heat
;****************************************

  q1 = dTdt - (omega*ss-advT)      ; term_1 - term_2
  q1@long_name = "q1: Net Heat Flux"
  q1@units     = "K/s"
  copy_VarCoords(T,q1)
  printVarSummary(q1) 
  printMinMax(q1,0)
  print("-----")

;****************************************
;---Apparent Heat Source
;****************************************

  Cp           = 1004.64 
  Cp@long_name = "specific heat of dry air at constant pressure"
  Cp@units     = "J/(K-kg)"  ; [kg-m2/s2]/(K-kg) => [kg-m2]/[K-kg-s2] => m2/[K-s2]  

  Q1  = Cp*q1  ; [J/(K-kg)][K/s] -> [J/(kg-s)] -> [(kg-m2/s2)/(kg-s)) -> [ m2/s3 ] 
               ; [J/(K-kg)][K/s] -> [(J/s)(1/kg)] -> W/kg   ???
  Q1@long_name = "Q1: Total Diabatic Heating as the Apparent Heat Source"
  Q1@units     = "m2/s3"     
  copy_VarCoords(T,Q1)
  printVarSummary(Q1) 
  printMinMax(Q1,0)
  print("-----") 

;*******************************************
;---Compute local dH/dt     
;*******************************************

  dHdt = center_finite_diff_n (H,time,False,0,0)
  copy_VarCoords(H, dHdt)
  dHdt@longname = "Specific Humidity: Local Time derivative"
  dHdt@units = "g/(kg-s)"    ; (g/kg)/s  
  printVarSummary(dHdt)
  printMinMax(dHdt,0)
  print("-----")

;****************************************
;---Compute advection term: global fixed grid: spherical harmonics
;---U*(dH/dlon) + V*(dH/dlat)
;****************************************

  long_name = "moisture advection"
  units     = "g/(kg-s)"     ; (m/s)*(g/kg)*(1/m)
  advH      = advect_variable(u,v,H,gridType,long_name,units,opt_adv) 

;****************************************
;---Compute vertical movement of H
;****************************************

  dHdp = center_finite_diff_n (H,p,False,1,npr)
  copy_VarCoords(H, dHdp)
  dHdp@longname = "Specific Humidity: Local Vertical Derivative"
  dHdp@units = "g/(kg-Pa)"   ; (g/kg)/Pa 
  printVarSummary(dHdp)
  printMinMax(dHdp,0)
  print("-----")

  dHdp  = omega*dHdp               ; overwrite ... convenience
  dHdp@longname = "Specific Humidity: Vertical Moisture Transport"
  dHdp@units    = "g/(kg-s)"       ; [(Pa/s)(g/kg)/Pa)]  
  printVarSummary(dHdp)
  printMinMax(dHdp,0)
  print("-----")
   
;****************************************
;---Apparent Moisture Sink
;****************************************

  q2    = -(dHdt +advH +dHdp)
  q2@long_name = "q2: moisture sink"      ; "?apparent? drying rate"
  q2@units     = "g/(kg-s)"
  copy_VarCoords(T,q2)
  printVarSummary(q2) 

  Lc           = 2.26e6    ; [J/kg]=[m2/s2]  Latent Heat of Condensation/Vaporization
  Lc@long_name = "Latent Heat of Condensation/Vaporization"
  Lc@units     = "J/kg"  ; ==> "m2/s2"

  Q2    = Lc*q2          ; (J/kg)(g/(kg-s))->(m2/s2)(g/(kg-s)) ->[(g/kg)][m2/s3]
                         ; J[oule]->  kg-m2/s2 = N-m = Pa/m3 = W-s       ; energy           
  Q2@long_name = "Q2: Apparent Moisture Sink"
  Q2@units     = "(g/kg)[m2/s3]"    
  copy_VarCoords(T,Q2)
  printVarSummary(Q2) 
  printMinMax(Q2,0)
  print("-----")
   
;****************************************
;---Make q1, q2 to per day
;****************************************

  q1 = q1*86400
  q2 = q2*86400
  q1@units     = "K/day"
  q2@units     = "g/(kg-day)"
   
;****************************************
;---Make Q1, Q2 to per W/m2
;****************************************

;;Q1 = Q1*?????
;;Q2 = Q2*?????
;;Q1@units = "W/m2"
;;Q2@units = "W/m2"

  return( [/q1, q2, Q1, Q2/] )
end
;==============================================================================
;                           MAIN
;==============================================================================
  netCDF  = True                                       ; Write netCDF 

;---Variable and file handling

  diri    = "/Users/shea/Data/netCDF/"
  f1      = addfile(diri+"temp_lag4tolead2.nc","r")    ; daily mean data
  f2      = addfile(diri+"omega_lag4tolead2.nc","r")     
  f3      = addfile(diri+"uwnd_lag4tolead2.nc","r")
  f4      = addfile(diri+"vwnd_lag4tolead2.nc","r")
  f5      = addfile(diri+"rhum_lag4tolead2.nc","r")
                                                       ; convenience: make all:  S->N
  temp    = f1->air(:,:,::-1,:)                        ; degK
  omega   = f2->omega(:,:,::-1,:)                      ; Pascal/s
  uwnd    = f3->uwnd(:,:,::-1,:)                       ; m/s
  vwnd    = f4->vwnd(:,:,::-1,:)                       ; m/s
  rhum    = short2flt(f5->rhum(:,:,::-1,:))            ; %   

;---Need specific humidity

  p       = f1->level                                  ; hPa [*]
  p4d     = conform(temp, p, 1)
  q       = mixhum_ptrh (p4d, temp, rhum,-2)           ; g/kg
  delete( [/p4d, rhum/] )
  q@long_name = "specific humidity"
  q@units = "g/kg"
  copy_VarCoords(temp, q)
  printVarSummary(q)
  printMinMax(q, 0)

;---Convert "hours since ..." to "seconds since ..."

  time    = f5->time                                   ; "hours since 1800-1-1"
  time    = time*3600                                          
  time@units  = "seconds since 1800-1-1 00:00:0.0"
  printVarSummary(time)
  print("---")
  t       = time                 ; Lyndz' name 

  ymdh = cd_calendar(time,-3)
  print(ymdh)

;---Convert hPa -> Pa for function

  p       = p*100                                      ; Pa  [100000,...,10000]
  p@units = "Pa"
  p!0     = "p"
  p&p     =  p                   ; not necessary
  printVarSummary(p)
  print("---")
  
;++++++++++++++++++++++++++++++++++++++++++++++++++++++
  Cp       = 1004.64      ; specific heat of dry air at constant pressure 
  Cp@units = "J/(K*kg)"  

  npr  = 1
  opt  = True
  q1q2 = Q1Q2_yanai(time,p,uwnd,vwnd,q,temp,omega,npr,opt)
  print(q1q2)
 
  q1   = q1q2[0]
  q2   = q1q2[1]
  printVarSummary(q1)
  printMinMax(q1,0)
  print("================")
  printVarSummary(q2)
  printMinMax(q2,0)
  print("================")

 ;Q1   = q1q2[2]
 ;Q2   = q1q2[3]

;********************************************
;---Vertical integration
; J[oule]      kg-m2/s2 = N-m = Pa/m3 = W-s       ; energy           
;********************************************

  ptop  = 30000.0                       ; Pa
  pbot  = 100000.0
  vopt  = 1 

  g     = 9.80665                       ; [m/s2] gravity at 45 deg lat used by the WMO
  dp    = dpres_plevel_Wrap(p,pbot,ptop,0)
  dpg   = dp/g                          ; Pa/(m/s2)=> (Pa-s2)/m   
  dpg@long_name = "Layer Mass Weighting"
  dpg@units     = "kg/m2"               ; dp/g     => Pa/(m/s2) => [kg/(m-s2)][m/s2] reduce to (kg/m2)
                                        ;             Pa (s2/m) => [kg/(m-s2)][s2/m]=>[kg/m2]
  dpg  := conform(q1,dpg,1)             ; reassign [convenience]

  q1int = wgt_vertical_n(q1,dpg,vopt,1) ; SUM[q1*dpg]  => (K/s)(kg/m2)
  q1int@long_name = "Heat Source: vertically integrated"
  q1int@units     = "(K/day)(kg/m2)"      ; (K/s)(kg/m2) => (kg-K)/(m2-s)
  printVarSummary(q1int)
  copy_VarCoords(temp(:,0,:,:),q1int)
  printMinMax(q1int,0)
  print("-----")

  q2int = wgt_vertical_n(q2,dpg,vopt,1)
  q2int@long_name = "Moisture Sink: vertically integrated"
  q2int@units     = "g/(day-m2)"
  copy_VarCoords(temp(:,0,:,:),q2int)
  printMinMax(q2int,0)
  print("-----")

;***********************************************
;---Save to a netcdf file
;***********************************************

  if (netCDF) 
      diro = "./"
      filo = "YANAI.apparent_heat.nc"
      ptho = diro+filo
      system("/bin/rm -f "+ptho)
      ncdf = addfile(ptho,"c")
    
      fAtt = True
      fAtt@title         = "Apparent Heat Source based on Yanai et al. 1973"
      fAtt@source_name   = "NCEP-NCAR Reanalysis 2"
      fAtt@source_URL    = "https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis2.html"
      fAtt@source        = "NOAA/OAR/ESRL PSD, Boulder, Colorado, USA"
      fAtt@Conventions   = "None"
      fAtt@creation_date = systemfunc("date")
      fileattdef(ncdf,fAtt)            ; copy file attributes
     
      filedimdef(ncdf,"time",-1,True) ; make time an UNLIMITED dimension
      ncdf->q1 = q1
      ncdf->q2 = q2
  end if

;***********************************************
;---Plot
;***********************************************

  nt    = 4
  YMDH  = ymdh(nt)
  LEVP  = 600

  opt = True
  opt@PrintStat = True
  stat_q1 = stat_dispersion(q1(nt,{LEVP},:,:), opt )
  stat_q2 = stat_dispersion(q2(nt,{LEVP},:,:), opt )
  stat_q1i= stat_dispersion(q1int(nt,:,:), opt )
  stat_q2i= stat_dispersion(q2int(nt,:,:), opt )

  plot  = new(2,graphic)

  wks   = gsn_open_wks("png","q2q1_yanai")        ; send graphics to PNG file
        
;--- mfc_adv and mfc_con at a specified pressure level

  res                   = True             ; plot mods desired
  res@gsnDraw           = False            ; don't draw yet
  res@gsnFrame          = False            ; don't advance frame yet

  res@cnFillOn          = True             ; turn on color
  res@cnLinesOn         = False            ; turn off contour lines
  res@cnLineLabelsOn    = False            ; turn off contour lines
 ;res@cnFillPalette     = "ViBlGrWhYeOrRe" ; set White-in-Middle color map
  res@cnFillPalette     = "amwg256"        ; set White-in-Middle color map
 ;res@cnFillMode        = "RasterFill"
  res@mpFillOn          = False            ; turn off map fill
;;res@lbLabelBarOn      = False            ; turn off individual cb's
  res@lbOrientation     = "Vertical"
                                           ; Use a common scale
  res@cnLevelSelectionMode = "ManualLevels"; manual set levels so lb consistent
  res@cnMaxLevelValF       =    5.0        ; max level
  res@cnMinLevelValF       = -res@cnMaxLevelValF     ; min level
  res@cnLevelSpacingF      =    0.20       ; contour interval

  res@gsnCenterString      = LEVP+"hPa"
  plot(0) = gsn_csm_contour_map(wks,q1(nt,{LEVP},:,:),res)

  res@cnMaxLevelValF       =    5.0        ; max level
  res@cnMinLevelValF       = -res@cnMaxLevelValF     ; min level
  res@cnLevelSpacingF      =    0.20       ; contour interval
  plot(1) = gsn_csm_contour_map(wks,q2(nt,{LEVP},:,:),res)

  resP                     = True          ; modify the panel plot
  resP@gsnMaximize         = True
  resP@gsnPanelMainString := YMDH
;;resP@gsnPanelLabelBar    = True          ; add common colorbar
  gsn_panel(wks,plot,(/2,1/),resP)         ; now draw as one plot

  delete(res@gsnCenterString) ; not used for this plot

  res@cnMaxLevelValF       =   14000.0     ; max level
  res@cnMinLevelValF       = -res@cnMaxLevelValF     ; min level
  res@cnLevelSpacingF      =    500.0      ; contour interval
  plot(0) = gsn_csm_contour_map(wks,q1int(nt,:,:),res)

  res@cnMaxLevelValF       =    7000.0     ; max level
  res@cnMinLevelValF       = -res@cnMaxLevelValF     ; min level
  res@cnLevelSpacingF      =    500.0      ; contour interval
  plot(1) = gsn_csm_contour_map(wks,q2int(nt,:,:),res)

  gsn_panel(wks,plot,(/2,1/),resP)         ; now draw as one plot


;---Cross section
 
  rescx                      = True                  ; plot mods desired
  rescx@gsnMaximize          = True

  LAT = 10 
  if (LAT.ge.0) then
      rescx@tiMainString     = YMDH+": "+LAT+"N"
  else
      rescx@tiMainString     = YMDH+": "+ABS(LAT)+"S"
  end if

  rescx@cnLevelSelectionMode = "ManualLevels"        ; manual contour levels
  rescx@cnLinesOn            = False            ; turn off contour lines
  rescx@cnLineLabelsOn       = False
  rescx@cnFillOn             = True                  ; turn on color fill
  rescx@cnFillPalette        = "ViBlGrWhYeOrRe" ; set White-in-Middle color map
  rescx@cnFillPalette        = "amwg256"        ; set White-in-Middle color map
  
  rescx@cnMaxLevelValF       =  5.0                  ; max level
  rescx@cnMinLevelValF       = -rescx@cnMaxLevelValF ; min level
  rescx@cnLevelSpacingF      =  0.25                 ; contour interval
  plt_q1 = gsn_csm_pres_hgt(wks,q1(nt,{1000:200},{LAT},:),rescx)
  rescx@cnMaxLevelValF       =  5.0                  ; max level
  rescx@cnMinLevelValF       = -rescx@cnMaxLevelValF ; min level
  rescx@cnLevelSpacingF      =  0.25                 ; contour interval
  plt_q2 = gsn_csm_pres_hgt(wks,q2(nt,{1000:200},{LAT},:),rescx)

具体Fortran代码

代码语言:javascript
复制
C*************************************************************** 
C     Compute voticity divegence             YANG 1992.10.25   
C   input UU  VV  OUTPUT VORT DIV (APPP)              * 
C*************************************************************** 
 program calsiq 
 parameter(I0=71,J0=71,K0=21,kk=19,k4=12,df1=1.,df2=1.,f0=-10.) 
      DIMENSION ww(i0,j0,kk,k4),w(i0,j0,kk),u(i0,j0,kk),v(i0,j0,kk) 
      DIMENSION uu(i0,j0,kk,k4),TT(i0,j0,kk,k4),HH(i0,j0,kk,k4),P(K0) 
     & ,t(i0,j0,kk),q(i0,j0,kk),QQ(i0,j0,kk,k4),vv(i0,j0,kk,k4)     
      DIMENSION PP(KK),siq1(i0,j0,kk,k4),siq2(i0,j0,kk,k4),h(i0,j0,kk) 
      DIMENSION K88(kk),work(i0,j0,k0,k4),OROG(i0,j0),hs(i0,j0,k4) 
      DIMENSION rq(i0,j0),tg(i0,j0,kk),th(i0,j0,kk),ft1(i0,j0,kk) 
      DIMENSION ft(i0,j0,kk),Q1(i0,j0),Q2(i0,j0),ft3(i0,j0,kk) 
c COMMON /CONTAA/RG,RD,CP,G,OMG,DLON,DLAT,FO,DF 
c     & ,DLATD,DLOND,CONV,COV,DPS12 
      CHARACTER*40,FILEN,FILEN1 
 DATA P/1000.,975.,950.,925.,900.,850.,800.,750.,700.,650.,600., 
     &  550.,500.,450.,400.,350.,300.,250.,200.,150.,100/  
 DATA PP/1000.,950.,900.,850.,800.,750.,700.,650.,600., 
     &  550.,500.,450.,400.,350.,300.,250.,200.,150.,100/  
c      DATA K66/1,2,3,4,5,9/ 
      DATA K88/1,3,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21/ 
C****************************************************************** 
C     Compute VORV DIAV    liuhz 1998.10.25           * 
C****************************************************************** 
 open(11,file='e:\data\ncep0629\data\orog1.grb' 
     & ,form='unformatted',access='direct',recl=i0*j0) 
  read(11,rec=1) 
     &((orog(I,J),I=1,i0),J=1,j0) 
 close(11) 
 open(11,file='e:\data\ncep0629\data\uu0629.grb' 
     & ,form='unformatted',access='direct',recl=i0*j0*k0*k4) 
  read(11,rec=1) 
     &((((work(I,J,k,K2),I=1,i0),J=1,j0),k=1,k0),k2=1,k4) 
 close(11) 
 do kt=1,k4 
 do k=1,kk 
 k2=k88(k) 
 do j=1,j0 
 do i=1,i0 
 uu(i,j,k,kt)=work(i,j,k2,kt) 
 enddo 
 enddo 
 enddo 
 enddo 
 open(11,file='e:\data\ncep0629\data\vv0629.grb' 
     & ,form='unformatted',access='direct',recl=i0*j0*k0*k4) 
  read(11,rec=1) 
     &((((work(I,J,k,K2),I=1,i0),J=1,j0),k=1,k0),k2=1,k4) 
 close(11) 
 do kt=1,k4 
 do k=1,kk 
 k2=k88(k) 
 do j=1,j0 
 do i=1,i0 
 vv(i,j,k,kt)=work(i,j,k2,kt) 
 enddo 
 enddo 
 enddo 
 enddo 
 open(11,file='e:\data\ncep0629\data\tt0629.grb' 
     & ,form='unformatted',access='direct',recl=i0*j0*k0*k4) 
  read(11,rec=1) 
     &((((work(I,J,k,K2),I=1,i0),J=1,j0),k=1,k0),k2=1,k4) 
 close(11) 
 do kt=1,k4 
 do k=1,kk 
 k2=k88(k) 
 do j=1,j0 
 do i=1,i0 
 tt(i,j,k,kt)=work(i,j,k2,kt) 
 enddo 
 enddo 
 enddo 
 enddo 
 open(11,file='e:\data\ncep0629\data\hh0629.grb' 
     & ,form='unformatted',access='direct',recl=i0*j0*k0*k4) 
  read(11,rec=1) 
     &((((work(I,J,k,K2),I=1,i0),J=1,j0),k=1,k0),k2=1,k4) 
 close(11) 
 do kt=1,k4 
 do k=1,kk 
 k2=k88(k) 
 do j=1,j0 
 do i=1,i0 
 hh(i,j,k,kt)=work(i,j,k2,kt) 
 enddo 
 enddo 
 enddo 
 enddo 
 open(11,file='e:\data\ncep0629\data\qq0629.grb' 
     & ,form='unformatted',access='direct',recl=i0*j0*k0*k4) 
  read(11,rec=1) 
     &((((work(I,J,k,K2),I=1,i0),J=1,j0),k=1,k0),k2=1,k4) 
 close(11) 
 do kt=1,k4 
 do k=1,kk 
 k2=k88(k) 
 do j=1,j0 
 do i=1,i0 
 qq(i,j,k,kt)=work(i,j,k2,kt) 
 enddo 
 enddo 
 enddo 
 enddo 
 irec=0 
 jrec=0 
 
 do 100 k=1,k4 
 do j=1,j0 
 do i=1,i0 
 do l=1,kk 
 h(i,j,l)=hh(i,j,l,k) 
 u(i,j,l)=uu(i,j,l,k) 
 v(i,j,l)=vv(i,j,l,k) 
 t(i,j,l)=tt(i,j,l,k) 
 q(i,j,l)=qq(i,j,l,k) 
 enddo 
 enddo 
 enddo 
 OPEN(20,FILE='e:\data\ncep0629\data\q1t.GRB',FORM='UNFORMATTED', 
     &ACCESS='DIRECT',RECL=i0*j0*(k0-1)) 
 OPEN(21,FILE='e:\data\ncep0629\data\q1v.GRB',FORM='UNFORMATTED', 
     &ACCESS='DIRECT',RECL=i0*j0*(k0-1)) 
 OPEN(22,FILE='e:\data\ncep0629\data\q1w.GRB',FORM='UNFORMATTED', 
     &ACCESS='DIRECT',RECL=i0*j0*(k0-1)) 
 OPEN(23,FILE='e:\data\ncep0629\data\w20.GRB',FORM='UNFORMATTED', 
     &ACCESS='DIRECT',RECL=i0*j0*(k0-1)) 
 OPEN(24,FILE='e:\data\ncep0629\data\q2t.GRB',FORM='UNFORMATTED', 
     &ACCESS='DIRECT',RECL=i0*j0*(k0-1)) 
 OPEN(25,FILE='e:\data\ncep0629\data\q2v.GRB',FORM='UNFORMATTED', 
     &ACCESS='DIRECT',RECL=i0*j0*(k0-1)) 
 OPEN(26,FILE='e:\data\ncep0629\data\q2w.GRB',FORM='UNFORMATTED', 
     &ACCESS='DIRECT',RECL=i0*j0*(k0-1)) 
 OPEN(27,FILE='e:\data\ncep0629\data\q12.GRB',FORM='UNFORMATTED', 
     &ACCESS='DIRECT',RECL=i0*j0) 
 
 DO L=1,kk 
 do i=1,i0 
 do j=1,j0 
 rq(i,j)=u(i,j,l) 
 enddo 
 enddo 
 call smth9(i0,j0,rq,rq) 
 do i=1,i0 
 do j=1,j0 
 u(i,j,l)=rq(i,j) 
 enddo 
 enddo 
 
 do i=1,i0 
 do j=1,j0 
 rq(i,j)=v(i,j,l) 
 enddo 
 enddo 
 call smth9(i0,j0,rq,rq) 
 do i=1,i0 
 do j=1,j0 
 v(i,j,l)=rq(i,j) 
 enddo 
 enddo 
 
 do i=1,i0 
 do j=1,j0 
 rq(i,j)=t(i,j,l) 
 enddo 
 enddo 
 call smth9(i0,j0,rq,rq) 
 do i=1,i0 
 do j=1,j0 
 t(i,j,l)=rq(i,j) 
 enddo 
 enddo 
 
 do i=1,i0 
 do j=1,j0 
 rq(i,j)=q(i,j,l) 
 enddo 
 enddo 
 call smth9(i0,j0,rq,rq) 
 do i=1,i0 
 do j=1,j0 
 q(i,j,l)=rq(i,j) 
 enddo 
 enddo 
 
 ENDDO 
 
 IF(K.GT.1)THEN 
 DO L=1,kk 
 do i=1,i0 
 do j=1,j0 
 rq(i,j)=tt(i,j,l,k-1) 
 enddo 
 enddo 
 call smth9(i0,j0,rq,rq) 
 do i=1,i0 
 do j=1,j0 
 tg(i,j,l)=rq(i,j) 
 enddo 
 enddo 
 do i=1,i0 
 do j=1,j0 
 rq(i,j)=qq(i,j,l,k-1) 
 enddo 
 enddo 
 call smth9(i0,j0,rq,rq) 
 do i=1,i0 
 do j=1,j0 
 ft1(i,j,l)=rq(i,j) 
 enddo 
 enddo 
 
      ENDDO 
 ENDIF 
 IF(K.LT.K4)THEN 
 DO L=1,kk 
 do i=1,i0 
 do j=1,j0 
 rq(i,j)=tt(i,j,l,K+1) 
 enddo 
 enddo 
 call smth9(i0,j0,rq,rq) 
 do i=1,i0 
 do j=1,j0 
 th(i,j,l)=rq(i,j) 
 enddo 
 enddo 
 
      do i=1,i0 
 do j=1,j0 
 rq(i,j)=qq(i,j,l,K+1) 
 enddo 
 enddo 
 call smth9(i0,j0,rq,rq) 
 do i=1,i0 
 do j=1,j0 
 ft3(i,j,l)=rq(i,j) 
 enddo 
 enddo 
 
 ENDDO 
 ENDIF 
 DT=24*3600.0             
  
C �����¶ȵľֵر仯 
 IF(K.EQ.1)THEN 
 DO I=1,i0 
 DO J=1,j0 
 DO L=1,kk 
 TG(I,J,L)=(TH(I,J,L)-T(I,J,L))/DT 
 fT(I,J,L)=(ft3(I,J,L)-q(I,J,L))/DT 
 ENDDO 
 ENDDO 
 ENDDO 
 ELSE IF(K.EQ.K4)THEN 
 DO I=1,i0 
 DO J=1,j0 
 DO L=1,kk 
 TG(I,J,L)=(T(I,J,L)-TG(I,J,L))/DT 
 fT(I,J,L)=(q(I,J,L)-fT1(I,J,L))/DT 
 ENDDO 
 ENDDO 
 ENDDO 
 ELSE 
 DO I=1,i0 
 DO J=1,j0 
 DO L=1,kk 
 TG(I,J,L)=(TH(I,J,L)-TG(I,J,L))/2.0/DT 
 fT(I,J,L)=(fT3(I,J,L)-fT1(I,J,L))/2.0/DT 
 ENDDO 
 ENDDO 
 ENDDO 
 ENDIF 
      call CALsiq12(U,V,t,q,h,tg,ft,OROG,irec,jrec 
     &          ,I0,j0,kk,df1,df2,f0,pp)  
 
 do l=1,kk 
 do j=1,j0 
 do i=1,i0 
c ampv(i,j,l,k)=ampv1(i,j,l,k)+ampv2(i,j,l,k) 
 enddo 
 enddo 
 enddo 
c      print*,ampv(12,18,11,kt) 
 WRITE(6,*)'TIME=',K 
 
c pause 1111 
100 continue 
 CLOSE(20);CLOSE(21);CLOSE(22);CLOSE(23);close(24);close(25) 
 close(26);close(27) 
 
cc      write(filen1,'(a)') 'e:\data\ncep0629\data\siq1.grb' 
ccc open(17,file=filen1, 
cc     &      form='unformatted',access='direct',recl=i0*j0*kk*k4) 
cc          write(17,rec=1) ((((siq1(I,J,K,kt),I=1,i0),J=1,j0) 
cc     &      ,K=1,kk),kt=1,k4) 
cc      CLOSE(17) 
cc      write(filen1,'(a)') 'e:\data\ncep0629\data\siq2.grb' 
cc open(17,file=filen1, 
cc     &      form='unformatted',access='direct',recl=i0*j0*kk*k4) 
cc          write(17,rec=1) ((((siq2(I,J,K,kt),I=1,i0),J=1,j0) 
cc     &      ,K=1,kk),kt=1,k4) 
cc      CLOSE(17) 
 
      STOP 
 END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
      SUBROUTINE CALsiq12(U,V,T,Q,h,tg,ft,OROG,irec,jrec 
     &                             ,I0,j0,k0,df1,df2,f0,pp)                 
C**************************************************************************      
C    Compute HEATIGH RATE (SOURCE/SINK) & WATER SOURCE/SINK    by liuhz 2009.3    *    
C    input  U=====> U WIND (m/s)                                                  *      
C           V=====> V WIND (m/s)                                                  * 
C           H ======> GEOPOTENTIAL HEIHGT (GEOPOTENTIAL METER) 
C           Q ====> SPECIFIC HUMIDITY (G/G)                                       * 
C           T =====> TEMPERATURE (K)                                              *      
C           tg =====> local variation of TEMPERATURE                              *      
C           ft =====> local variation of SPECIFIC HUMIDITY  
C         OROG======> OROGRAPHIC HEIHGT(GEOPOTENTIAL METER)                       *      
C    output UG(Q1V)  =====> ADVECTIVE TRANSFER of HEAT (K/S)  
C           WG(Q1W)  =====> VERTICAL TRANSFER  OF HEAT (K/S)      
C           W(W20)   =====> W VELOCITY (HPA/S) ;                                  *      
C           FV(Q2V)  =====> ADVECTIVE TRANSFER of water (G/G/S)  
C           FW(Q2W)  =====> VERTICAL TRANSFER  OF water (G/G/S)      
C           Q1(Q12)  =====> HEAT SOURCE/SINK (K/S) ;       
C           Q2       =====> WATER SOURCE/SINK (G/G/S) ;   
C           tg (Q1T) =====> local variation of TEMPERATURE (K/S)                  *      
C           ft (Q2T) =====> local variation of SPECIFIC HUMIDITY (G/G/S) 
C**************************************************************************      
C      INCLUDE 'PARAMET213' 
C      PARAMETER (I0=161,J0=113,K0=13,DF1=1.0,DF2=1.0,F0=-10.)           
C       I0,J0,K0---GRID NUMBER IN X Y Z; DF1-DF2- FIELD RESOLUTION IN X Y;       
C       F0--FIRST LATITUDE FROM SOUTH TO NORTH                                   
C       P(1) --P(K0) EACH LAYEL PRESSURE FROM LOW TO HIHG                        
C      DATA PL/1000.0,925.0,850.0,700.0,600.,500.0,400.0,300.0,250.0,200.        
C     &   ,150.,100.,50./                                                        
C ���ò�ֵ��NCEP/NCAR 20��U��V��T��H���ϼ���Q1 & Q2 
C ����X��Y��P��T�������� 
 dimension OROG(i0,j0),H(i0,j0,k0),PP(k0),a(i0,j0,k0) 
     &          ,siq1(I0,J0,K0),siq2(I0,J0,K0),AR(I0,J0,K0),P(k0)                         
      DIMENSION U(I0,J0,K0),V(I0,J0,K0),PL(K0),T(I0,J0,K0)                 
 DIMENSION PS(i0,j0),US(i0,j0),VS(i0,j0),WS(i0,j0),ts(i0,j0) 
 DIMENSION ILAY(i0,j0),q(i0,j0,k0),gtmp(i0,j0,k0) 
 DIMENSION DIVG(i0,j0,k0),DIVS(i0,j0),W(i0,j0,k0) 
 dimension wg(i0,j0,k0),ug(i0,j0,k0),tg(i0,j0,k0),TH(i0,j0,k0) 
 DIMENSION rQ(i0,j0),HT(i0,j0,k0),wt(i0,j0),Q1(i0,j0),Q2(i0,j0) 
 dimension ft(i0,j0,k0),fv(i0,j0,k0),fw(i0,j0,k0) 
 dimension ft1(i0,j0,k0),ft3(i0,j0,k0),ff(i0,j0,k0) 
 
      COV=3.1415926/180.                                                         
c      DF=DF2                                                                    
      DLON=DF1*COV                                                               
      DLAT=DF2*COV                                                               
      DLOND=DF1                                                                  
      DLATD=DF2                                                                  
      RG=6371229.0                                                               
      DY=2.0*RG*DLAT   
 ck=0.286 
 gama=0.0065 
 rd=287.04 
 g=9.80616 
 dk=g/(rd*gama) 
 DP=-50.0 
 WTt=0.0 
c DT=6*3600.0             
 DT=24*3600.0             
 CP=1005.0 
 LH=2501000.0 
  
C �����¶�ƽ��---������ 
 do i=2,i0-1 
 do j=2,j0-1 
      F1=(F0+FLOAT(J-1)*DLATD)*COV                                               
      DX=2.0*RG*COS(F1)*DLON                                                     
  DO L=1,k0 
  UG(I,J,L)=U(I,J,L)*(T(I+1,J,L)-T(I-1,J,L))/DX 
c     &       COS((-87.0+1.0*(j-1))*RAD)+ 
     &     +V(I,J,L)*(T(I,J+1,L)-T(I,J-1,L))/DY 
  fv(I,J,L)=U(I,J,L)*(q(I+1,J,L)-q(I-1,J,L))/DX 
c     &       COS((-87.0+1.0*(j-1))*RAD)+ 
     &     +V(I,J,L)*(q(I,J+1,L)-q(I,J-1,L))/DY 
  ENDDO 
 ENDDO 
 ENDDO 
 
 do j=2,j0-1 
      F1=(F0+FLOAT(J-1)*DLATD)*COV                                               
      DX=2.0*RG*COS(F1)*DLON                                                     
  DO L=1,k0 
  UG(1,J,L)=U(1,J,L)*(T(2,J,L)-T(i0-1,J,L))/DX 
c     &       COS((-87.0+1.0*(j-1))*RAD)+ 
     &     +V(1,J,L)*(T(1,J+1,L)-T(1,J-1,L))/DY 
  fv(1,J,L)=U(1,J,L)*(q(2,J,L)-q(i0-1,J,L))/DX 
c     &       COS((-87.0+1.0*(j-1))*RAD)+ 
     &     +V(1,J,L)*(q(1,J+1,L)-q(1,J-1,L))/DY 
  ENDDO 
  DO L=1,k0 
  UG(i0,J,L)=U(i0,J,L)*(T(1,J,l)-T(i0-1,J,L))/DX 
c     &       COS((-87.0+1.0*(j-1))*RAD)+ 
     &     +V(i0,J,L)*(T(i0,J+1,L)-T(i0,J-1,L))/DY 
  fv(i0,J,L)=U(i0,J,L)*(q(1,J,l)-q(i0-1,J,L))/DX 
c     &       COS((-87.0+1.0*(j-1))*RAD)+ 
     &     +V(i0,J,L)*(q(i0,J+1,L)-q(i0,J-1,L))/DY 
  enddo 
 ENDDO 
 do l=1,k0 
  CALL BOU1(UG(1,1,L),I0,J0) 
  CALL BOU1(FV(1,1,L),I0,J0) 
 enddo 
c      pause 888 
  
 DO I=1,i0 
 DO J=1,j0 
C ����� 
  DO L=1,k0 
  gtmp(i,j,l)=t(i,j,l)*(1000.0/pp(l))**ck 
  ENDDO 
c ���㶥�㴹ֱ�ٶ� 
  WT(I,J)=-(TG(I,J,k0)+UG(I,J,k0))/ 
     &  ((GTMP(I,J,k0)-GTMP(I,J,k0-1))/DP*(PP(k0)/1000.0)**CK) 
 
C ������ڵ��ε���Ͳ� 
c  orog(i,j)=orog(i,j)*g 
  if(h(i,j,1).gt.orog(i,j))then 
  ip=1 
  else 
  DO L=1,k0-1 
  IF((H(I,J,L).LE.OROG(I,J)).AND. 
     & (H(I,J,L+1).GT.OROG(I,J)))THEN 
   IP=L+1 
  ELSE 
  ENDIF 
  ENDDO 
  endif 
   ILAY(I,J)=IP 
       
C ������¡�������ѹ---��Ԫ������� 
  ts(i,j)=t(i,j,ip)-gama*(orog(i,j)-h(i,j,ip)) 
C �ڲ������ѹ������� 
  if(ip.gt.1)then 
  CC=(H(I,J,IP)-OROG(I,J))/(H(I,J,IP)-H(I,J,IP-1)) 
  CD=(PP(IP-1)/PP(IP))**CC 
  PS(I,J)=CD*PP(IP) 
  US(I,J)=U(I,J,IP)+CC*(U(I,J,IP)-U(I,J,IP+1)) 
  VS(I,J)=V(I,J,IP)+CC*(V(I,J,IP)-V(I,J,IP+1)) 
  else 
  ps(i,j)=pp(ip)*(1-gama*(orog(i,j)-h(i,j,ip))/t(i,j,ip))**dk 
  us(i,j)=u(i,j,ip)+(U(I,J,IP)-U(I,J,IP+1))* 
     &  (OROG(I,J)-H(I,J,IP))/(H(I,J,IP)-H(I,J,IP+1)) 
  US(I,J)=0.6*US(I,J) 
  Vs(i,j)=V(i,j,ip)+(V(I,J,IP)-V(I,J,IP+1))* 
     &  (OROG(I,J)-H(I,J,IP))/(H(I,J,IP)-H(I,J,IP+1)) 
  VS(I,J)=0.6*VS(I,J) 
  endif 
 enddo 
 enddo 
C ������洹ֱ�ٶ� 
 do i=2,i0-1 
 do j=2,j0-1 
      F1=(F0+FLOAT(J-1)*DLATD)*COV                                               
      DX=2.0*RG*COS(F1)*DLON                                                     
 
  WS(I,J)=(-1.0)*PS(I,J)*g/Rd/Ts(I,J)* 
     &   (US(I,J)*(OROG(I+1,J)-OROG(I-1,J))/DX 
c     &       COS((-87.0+1.0*(j-1))*RAD)+ 
     &     +VS(I,J)*(OROG(I,J+1)-OROG(I,J-1))/DY) 
 ENDDO 
 ENDDO 
  CALL BOU1(WS,I0,J0) 
 
c do j=2,j0-1 
c  WS(1,J)=(-1.0)*PS(1,J)*g/Rd/Ts(1,J)* 
c     &   US(1,J)*(OROG(2,J)-OROG(i0-1,J))/DX 
ccc     &       COS((-87.0+1.0*(j-1))*RAD)+ 
c     &     +VS(I,J)*(OROG(I,J+1)-OROG(I,J-1))/DY 
c  WS(i0,J)=(-1.0)*PS(i0,J)*g/Rd/Ts(i0,J)* 
c     &   (US(i0,J)*(OROG(1,J)-OROG(i0-1,J))/DX 
c     &       COS((-87.0+1.0*(j-1))*RAD)+ 
c     &     +VS(I,J)*(OROG(I,J+1)-OROG(I,J-1))/DY) 
c enddo 
C �������ɢ�� 
 DO I=2,i0-1 
 DO J=2,j0-1 
      F1=(F0+FLOAT(J-1)*DLATD)*COV                                               
      DX=2.0*RG*COS(F1)*DLON                                                     
      TGA=TAN(F1)/RG 
 
 DO L=1,k0 
  DIVG(I,J,L)=(U(I+1,J,L)-U(I-1,J,L))/DX 
c     &       /COS((-87.0+1.0*(j-1))*RAD)/(2*DD)+ 
     &    +(V(I,J+1,L)-V(I,J-1,L))/DY-V(I,J,L)*TGA 
 ENDDO 
 ENDDO 
 ENDDO 
 iI=1 
 DO J=2,j0-1 
      F1=(F0+FLOAT(J-1)*DLATD)*COV                                               
      DX=2.0*RG*COS(F1)*DLON                                                     
      TGA=TAN(F1)/RG 
 
 DO L=1,k0 
  DIVG(II,J,L)=(U(II+1,J,L)-U(i0,J,L))/DX 
     &   +(V(II,J+1,L)-V(II,J-1,L))/DY-V(II,J,L)*TGA 
 ENDDO 
 ENDDO 
 iI=i0 
 DO J=2,j0-1 
      F1=(F0+FLOAT(J-1)*DLATD)*COV                                               
      DX=2.0*RG*COS(F1)*DLON                                                     
      TGA=TAN(F1)/RG 
 
 DO L=1,k0 
 DIVG(II,J,L)=(U(1,J,L)-U(II-1,J,L))/DX 
c     &       /COS((-87.0+1.0*(j-1))*RAD)/(2*DD)+ 
     &   +(V(II,J+1,L)-V(II,J-1,L))/DY-V(II,J,L)*TGA 
 ENDDO 
 ENDDO 
 do l=1,k0 
  CALL BOU1(DIVG(1,1,l),I0,J0) 
 enddo 
 
C �������ɢ�� 
 DO I=2,i0-1 
 DO J=2,j0-1 
      F1=(F0+FLOAT(J-1)*DLATD)*COV                                               
      DX=2.0*RG*COS(F1)*DLON                                                     
      TGA=TAN(F1)/RG 
 
  DIVS(I,J)=(US(I+1,J)-US(I-1,J))/DX 
     &  +(VS(I,J+1)-VS(I,J-1))/DY-VS(I,J)*TGA 
 ENDDO 
 ENDDO 
 iI=1 
 DO J=2,j0-1 
      F1=(F0+FLOAT(J-1)*DLATD)*COV                                               
      DX=2.0*RG*COS(F1)*DLON                                                     
      TGA=TAN(F1)/RG 
 
  DIVS(II,J)=(US(II+1,J)-US(i0,J))/DX 
     &  +(VS(II,J+1)-VS(II,J-1))/DY-VS(II,J)*TGA 
 ENDDO 
 iI=i0 
 DO J=2,j0-1 
      F1=(F0+FLOAT(J-1)*DLATD)*COV                                               
      DX=2.0*RG*COS(F1)*DLON                                                     
      TGA=TAN(F1)/RG 
 
  DIVS(II,J)=(US(1,J)-US(II-1,J))/DX 
c     &       /COS((-87.0+1.0*(j-1))*RAD)/(2*DD)+ 
     &   +(VS(II,J+1)-VS(II,J-1))/DY-VS(II,J)*TGA 
 ENDDO 
  CALL BOU1(DIVS,I0,J0) 
 
C ���㴹ֱ�ٶ� 
 do i=1,i0 
 do j=2,j0-1 
  IP=ILAY(I,J) 
   if(ip.GT.1)THEN 
   DO L=1,IP-1 
   W(I,J,L)=-9.99E33 
   ENDDO 
   ELSE 
   ENDIF 
  W(I,J,IP)=WS(I,J)-(DIVS(I,J)+DIVG(I,J,IP))/2*(PP(IP)-PS(I,J)) 
  do L=IP+1,k0 
   W(I,J,L)=W(I,J,L-1)-(DIVG(I,J,L-1)+DIVG(I,J,L))/2*DP 
  ENDDO 
   MM=(k0-IP+1)*(k0-IP+2) 
  DO L=IP,k0 
   LL=L-IP+1 
   W(I,J,L)=W(I,J,L)-LL*(LL+1)*(W(I,J,k0)-WT(i,j))/REAL(MM) 
  ENDDO 
C ����λ�´�ֱ���� 
C  if(ip.gt.1)then 
C  do l=1,ip-1 
C  wg(i,j,l)=-9.99e33 
C  enddo 
C  else 
C  endif 
  do l=ip+1,k0-1 
  wg(i,j,l)=W(I,J,L)*(GTMP(I,J,L+1)-GTMP(I,J,L-1))/2/DP 
  fw(i,j,l)=W(I,J,L)*(q(I,J,L+1)-q(I,J,L-1))/2/DP 
  ENDDO 
  WG(I,J,IP)=W(I,J,IP)*(GTMP(I,J,IP+1)-GTMP(I,J,IP))/DP 
  WG(I,J,k0)=W(I,J,k0)*(GTMP(I,J,k0)-GTMP(I,J,k0-1))/DP 
  fW(I,J,IP)=W(I,J,IP)*(q(I,J,IP+1)-q(I,J,IP))/DP 
  fW(I,J,k0)=W(I,J,k0)*(q(I,J,k0)-q(I,J,k0-1))/DP 
 
 ENDDO 
 ENDDO 
 do L=ip+1,k0-1 
  CALL BOU1(WG(1,1,L),I0,J0) 
  CALL BOU1(fW(1,1,L),I0,J0) 
      enddo 
 
C ��������� 
 DO I=1,i0 
 DO J=2,j0-1 
  IP=ILAY(I,J) 
  IF(IP.GT.1)THEN 
  DO L=1,IP-1 
   HT(I,J,L)=9.99e20 
   ff(i,j,l)=9.99e20 
   tg(i,j,l)=9.99e20 
   ug(i,j,l)=9.99e20 
   wg(i,j,l)=9.99e20 
   ft(i,j,l)=9.99e20 
   fv(i,j,l)=9.99e20 
   fw(i,j,l)=9.99e20 
  ENDDO 
  ELSE 
  ENDIF 
   DO L=IP,k0 
  wg(i,j,l)=WG(I,J,L)*(PP(L)/1000.0)**CK 
  ft(i,j,l)=-ft(i,j,l)*lh/cp 
  fv(i,j,l)=-fv(i,j,l)*lh/cp 
  fw(i,j,l)=-fw(i,j,l)*lh/cp 
   HT(I,J,L)=TG(I,J,L)+UG(I,J,L)+WG(I,J,L) 
   ff(I,J,L)=fT(I,J,L)+fv(I,J,L)+fW(I,J,L) 
   ENDDO 
  
 SS=0.0 
 sx=0.0 
 do L=IP,k0-1 
 SI=(HT(I,J,L)+HT(I,J,L+1))*DP/2.0 
 SS=SS-SI 
 Sj=(ff(I,J,L)+ff(I,J,L+1))*DP/2.0 
 Sx=Sx-Sj 
 ENDDO 
 
 Q1(I,J)=SS/G*100.0*CP 
 q2(i,j)=sx/g*100.0*cp 
 
 ENDDO 
 ENDDO 
 do L=ip,k0-1 
  CALL BOU1(HT(1,1,L),I0,J0) 
  CALL BOU1(ff(1,1,L),I0,J0) 
      enddo 
  CALL BOU1(q1,I0,J0) 
  CALL BOU1(q2,I0,J0) 
 
 irec=irec+1 
 write(20,REC=irec)(((tg(I,J,L),I=1,i0),J=1,j0),l=1,k0-1) 
 write(21,REC=irec)(((ug(I,J,L),I=1,i0),J=1,j0),l=1,k0-1) 
 write(22,REC=irec)(((wg(I,J,L),I=1,i0),J=1,j0),l=1,k0-1) 
 write(23,REC=irec)(((w(I,J,L),I=1,i0),J=1,j0),l=1,k0-1) 
 write(24,REC=irec)(((ft(I,J,L),I=1,i0),J=1,j0),l=1,k0-1) 
 write(25,REC=irec)(((fv(I,J,L),I=1,i0),J=1,j0),l=1,k0-1) 
 write(26,REC=irec)(((fw(I,J,L),I=1,i0),J=1,j0),l=1,k0-1) 
 
 jrec=jrec+1 
 write(27,REC=jrec)((q1(I,J),I=1,i0),J=1,j0) 
 jrec=jrec+1 
 write(27,REC=jrec)((q2(I,J),I=1,i0),J=1,j0) 
 
      print*,irec,jrec 
1999  continue 
 
 return 
 end 
ccccccccccccccccccccccccccccccccccccccccccccccccccccc 
 subroutine smth9(i0,j0,ff,gg) 
 real ff(i0,j0),gg(i0,j0) 
 do i=2,i0-1 
 do j=2,j0-1 
 s=0.5 
 gg(i,j)=ff(i,j)+s/2.0*(1-s)*(ff(i+1,j)+ff(i-1,j)+ff(i,j+1) 
     & +ff(i,j-1)-4*ff(i,j))+s*s/4*(ff(i+1,J+1)+ff(i-1,j-1)+ 
     &  ff(i+1,j-1)+ff(i-1,j+1)-4*ff(i,j)) 
 enddo 
 enddo 
 i=1 
 do j=2,j0-1 
 gg(i,j)=ff(i,j)+s/2.0*(1-s)*(ff(i+1,j)+ff(i0,j)+ff(i,j+1) 
     & +ff(i,j-1)-4*ff(i,j))+s*s/4*(ff(i+1,J+1)+ff(i0,j-1)+ 
     &  ff(i+1,j-1)+ff(i0,j+1)-4*ff(i,j)) 
 enddo 
 i=i0 
 do j=2,j0-1 
 gg(i,j)=ff(i,j)+s/2.0*(1-s)*(ff(1,j)+ff(i-1,j)+ff(i,j+1) 
     & +ff(i,j-1)-4*ff(i,j))+s*s/4*(ff(1,J+1)+ff(i-1,j-1)+ 
     &  ff(i,j-1)+ff(i-1,j+1)-4*ff(i,j)) 
 enddo 
 j=1 
 do i=1,i0 
 gg(i,j)=ff(i,j) 
 enddo 
 j=j0 
 do i=1,i0 
 gg(i,j)=ff(i,j) 
 enddo 
 return 
 end 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
      SUBROUTINE BOU1(AA,MX,MY) 
      DIMENSION AA(MX,MY) 
      SS=.125 
      M=0 
      DO 1 I=2,MX-1 
      AA(I,1)=AA(I,2) 
      AA(I,MY)=AA(I,MY-1) 
 1    CONTINUE 
      DO 2 J=2,MY-1 
      AA(1,J)=AA(2,J) 
      AA(MX,J)=AA(MX-1,J) 
 2    CONTINUE 
      AA(1,1)=AA(2,1) 
      AA(1,MY)=AA(2,MY) 
      AA(MX,1)=AA(MX,2) 
      AA(MX,MY)=AA(MX,MY-1) 
      RETURN 
      END 

IDL代码

Q1Q2_Y73: Calculate apparent heat source and moisture sink, after Yanai et al. (1973). Files required:

代码语言:javascript
复制
div_vel_sph.pro
omega_vel.pro
press2alt.pro (by Dominik Brunner)
q1q2_y73.pro

地址:http://www.johnny-lin.com/lib.shtml

参考:

  1. http://www.ncl.ucar.edu/Applications/wind.shtml2.
  2. http://read.pudn.com/downloads378/sourcecode/others/1631672/calswa.f__.htm
  3. http://www.johnny-lin.com/lib.shtml
本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2020-12-14,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 气象杂货铺 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体分享计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 具体NCL代码如下:
  • IDL代码
相关产品与服务
对象存储
对象存储(Cloud Object Storage,COS)是由腾讯云推出的无目录层次结构、无数据格式限制,可容纳海量数据且支持 HTTP/HTTPS 协议访问的分布式存储服务。腾讯云 COS 的存储桶空间无容量上限,无需分区管理,适用于 CDN 数据分发、数据万象处理或大数据计算与分析的数据湖等多种场景。
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档