前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >气象编程 | 水汽通量:水汽通量散度(Moisture Flux Convergence)以及垂直积分

气象编程 | 水汽通量:水汽通量散度(Moisture Flux Convergence)以及垂直积分

作者头像
MeteoAI
发布2020-12-22 15:19:03
11.2K0
发布2020-12-22 15:19:03
举报
文章被收录于专栏:MeteoAIMeteoAI

mfc_div_1.ncl: Calculate various divergence and moisture quantities including Vertically Integrated Moisture Flux Convergence (VIMFC). VIMFC has a high correlation with frontal and convective activity. Positive values indicate net precipitation. The following equation is implemented within mfc_div_1.ncl

This example uses uv2dvF_Wrap [uv2dvF] because the grid is a global fixed grid. For global gaussian, uv2dvG_Wrap [uv2dvG] should be used. For a regional grid uv2dv_cfd should be used.

mfc_div_1.ncl

代码语言:javascript
复制
;*************************************************
; mfc_div_1.ncl
;
; Concepts illustrated:
;   - Read daily mean wind components, humidity and sfc. pressure 
;     from different files
;   - Reorder the input (N==>S) grid order to (S==>N) via NCL syntax  ::-1
;   - Calculate mass weightined layer thickness [units="kg/m2"]
;   - Calculate moisture flux [uq, vq]             
;   - Calculate moisture flux divergence using spherical harmonics              
;   - Integrate the moisture flux divergence using mass weighting 
;   - Plot a number of quantities
;*************************************************
;---Calculate the Horizontal Moisture Flux Convergence [MFC]
;*************************************************
;---High frequency source data: hourly/3hr/6hr/12hr/daily .... NOT monthly values
;---References:
;---http://www.cgd.ucar.edu/cas/catalog/newbudgets/
;---http://tornado.sfsu.edu/geosciences/classes/e260/AtmosphericRivers/Moisture%20Flux.pdf
;---https://www.spc.noaa.gov/publications/banacos/mfc-sls.pdf
;===================================================================
;   Data Source: ESRL Physical Sciences Division
;        https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html
;   NCEP Reanalysis data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, 
;   from their Web site at https://www.esrl.noaa.gov/psd/
;===================================================================

  ptop = 300              ; 'shum' upper level                  
  ptop@units = "hPa"
  g    = 9.80665          ; m/s2

  date  = 20080715        ; NH summer

;---ESRL: CDC data

  diri = "./"
  filq = "shum.2008.nc"   ; daily data for current year [366 days]
  filu = "uwnd.2008.nc"
  filv = "vwnd.2008.nc"
  filps= "pres.sfc.2008.nc"

  pthu = diri+filu      
  pthv = diri+filv
  pthq = diri+filq     
  pthps= diri+filps   

  fu   = addfile(pthu ,"r")
  fv   = addfile(pthv ,"r")
  fq   = addfile(pthq ,"r")
  fps  = addfile(pthps,"r")

;---Time

  ymd  = cd_calendar(fu->time, -2)    ; ymd[*]: human readable
  nt   = ind(ymd.eq.date)             ; date for plotting and testing
 
TEST = True        
if (.not.TEST) then                   ; all times                   
  u    = fu->uwnd(:,{1000:ptop},:,:)  ; m/s, (time,level,lat,lon)
  v    = fv->vwnd(:,{1000:ptop},:,:)
  q    = fq->shum                     ; [kg/kg], 1000-300 levels only
  ps   = fps->pres                    ; Pa=>[kg/(m-s2)], (time,lat,lon)

else                                  ; one time step; keep time dimension [ nt:nt: ]                     
  u    = fu->uwnd(nt:nt,{1000:ptop},:,:); m/s, (time,level,lat,lon)
  v    = fv->vwnd(nt:nt,{1000:ptop},:,:)
  q    = fq->shum(nt:nt,:,:,:)        ; [kg/kg], 1000-300 levels only
  ps   = fps->pres(nt:nt,:,:)         ; Pa=>[kg/(m-s2)], (time,lat,lon)
  nt   = 0                            ; only one time step
end if

;---Vertical levels

  ptop = ptop*100
  ptop@units = "Pa"

  plev = q&level                      ; hPa  
  plev = plev*100                     ; [100000,...,30000] Pa [kg/(m-s2)]
  plev@units = "Pa"

;---Change [kg/kg] to [g/kg]; not necessary: but common units for q

  q    = q*1000            
  q@units = "g/kg"

;---Divergence function [used later] requires S->N grid order

  u    = u(:,:,::-1,:)  
  v    = v(:,:,::-1,:)
  q    = q(:,:,::-1,:)     
  ps   =ps(:,  ::-1,:)       

;---Layer thickness: ; Pa=>[kg/(m-s2)], (time,level,lat,lon) 
;---Mass weighting:  (dp/g) => [Pa/(m/s2)] => (Pa-s2)/m => [kg/(m-s2)][s2/m] =>  (kg/m2)
;---Reference: http://www.cgd.ucar.edu/cas/catalog/newbudgets/

  dp   = dpres_plevel_Wrap(plev, ps, ptop, 0) ; Pa; layar thickness 

  dpg  = dp/g    
  dpg@long_name = "Layer Mass Weighting"
  dpg@units     = "kg/m2"                     ; dp/g, Pa/(m s-2), reduce to kg m-2

;---Moisture flux components at each pressure level

  uq   = u*q                                  ; (:,:,:,:)
  uq@long_name = "Zonal Moisture Flux [uq]"
  uq@units = "["+u@units+"]["+q@units+"]"     ; [m/s][g/kg]     
  copy_VarCoords(u,uq)                        ; (time,level,lat,lon)

  vq   = v*q                                  ; (:,:,:,:)
  vq@long_name = "Meridional Moisture Flux [vq]"
  vq@units = "["+v@units+"]["+q@units+"]" 
  copy_VarCoords(v,vq)                        ; (time,level,lat,lon)

PRINT_RAW = True
if (PRINT_RAW) then
  printVarSummary(q)                          ; (time,level,lat,lon); g/kg
  printMinMax(q,0)
  print("-----")
  printVarSummary(u)                          ; (time,level,lat,lon); m/s
  printMinMax(u,0)
  print("-----")
  printVarSummary(v)
  printMinMax(v,0)
  print("-----")
  printVarSummary(ps)                         ; (time,lat,lon); Pa => kg/(m-s2) 
  printMinMax(ps,0)
  print("-----")
  printVarSummary(uq)                         ; (time,level,lat,lon); (m/s)(g/kg)
  printMinMax(uq,0)
  print("-----")
  printVarSummary(vq)
  printMinMax(vq,0)
  print("-----")
  printVarSummary(dp)                         ; (time,level,lat,lon); Pa => kg/(m-s2)
  printMinMax(dp,0)
  print("-----")
                                ; examine layer thickness at selected locations
  print(dp(nt,:,{40},{180}))    ; mid-Pacific
  print(dp(nt,:,{40},{255}))    ; Boulder, CO 
  print("-----")
end if

;---Integrated mass weighted moisture flux components

  uq_dpg = uq*dpg                ; mass weighted 'uq'; [m/s][g/kg][kg/m2]=>[m/s][g/kg]
  iuq    = dim_sum_n(uq_dpg, 1)
  iuq@long_name = "Integrated Zonal UQ [uq*dpg]" 
  iuq@LONG_NAME = "Sum: Mass Weighted Integrated Zonal Moisture Flux [uq*dpg]" 
  iuq@units     = "[m/s][g/kg]"
  copy_VarCoords(u(:,0,:,:), iuq); (time,lat,lon)
  delete(uq_dpg)

  vq_dpg = vq*dpg                ; mass weighted 'vq'; [m/s][g/kg][kg/m2]=>[m/s][g/kg] 
  ivq    = dim_sum_n(vq_dpg, 1)
  ivq@long_name = "Integrated Meridional VQ [vq*dpg]" 
  ivq@LONG_NAME = "Sum: Mass Weighted Integrated Meridional Moisture Flux [vq*dpg]" 
  ivq@units     = "[m/s][g/kg]"
  copy_VarCoords(v(:,0,:,:), ivq); (time,lat,lon)
  delete(vq_dpg)

;---Divergence of moisture flux: uv2dvF => global 'fixed' rectilinear grid

  duvq  = uv2dvF_Wrap(uq, vq)    ; (time,level,lat,lon)
  duvq@long_name = "Divergence of Moisture Flux"
  duvq@units     = "g/(kg-s)"    ; (1/m)*[(m/s)(g/kg)] => [g/(kg-s)]

;---Mass weighted integration [sum] of the divergence of moisture flux

  duvq_dpg = duvq*dpg            ;  [g/(kg-s)][kg/m2] => [g/(m2-s)]
  iduvq    = dim_sum_n(duvq_dpg, 1)
  iduvq@long_name = "Integrated Mass Wgt MFC" 
  iduvq@LONG_NAME = "Integrated Mass Weighted Moisture Flux Convergence" 
  iduvq@units     = "g/(m2-s)"
  copy_VarCoords(u(:,0,:,:), iduvq)      ; (time,lat,lon)
  delete(duvq_dpg)

  VIMFC =  iduvq           ; keep meta data                         
  VIMFC = -VIMFC           ; Note the preceding -1 [negative precedes integration] 
  VIMFC@long_name = "VIMFC"

;---Another way to compute Integrated divergence of moisture flux [iduvq_1]

;;IDUVQ = wgt_vertical_n(duvq, dp, 2, 1) 
;;iduvq_0 = IDUVQ[0]
;;iduvq_0 = iduvq_0/g                    ; complete mass weighting
;;iduvq_0@long_name = "Average Mass Weighted MFC" 
;;iduvq_0@LONG_NAME = "Average Mass Weighted Moisture Flux Convergence" 
;;iduvq_0@units     = "g/(m2-s)"

;;iduvq_1 = IDUVQ[1]                     ; same as iuvq_sum
;;iduvq_1 = iduvq_1/g
;;iduvq_1@long_name = "Integrated MFC" 
;;iduvq_1@LONG_NAME = "Integrated Moisture Flux Convergence" 
;;iduvq_1@units     = "g/(m2-s)"

PRINT_RESULT = True
if (PRINT_RESULT) then
    printVarSummary(iuq)                 ; (time,lat,lon)
    printMinMax(iuq,0)
    print("-----")
    printVarSummary(ivq)                 ; (time,lat,lon)
    printMinMax(ivq,0)
    print("-----")
    printVarSummary(duvq)                ; (time,lev,lat,lon)
    printMinMax(duvq,0)
    print("-----")
    printVarSummary(iduvq)               ; (time,lat,lon)
    printMinMax(iduvq,0)
    print("-----")
;;  printVarSummary(iduvq_0)             ; (time,lat,lon)
;;  printMinMax(iduvq_0,0)
;;  print("-----")
;;  printVarSummary(iduvq_1)             ; (time,lat,lon)
;;  printMinMax(iduvq_1,0)
;;  print("-----")
end if

;*************************************************
; Calculate divergence: Use Wrap to include meta data
; Calculate divergent wind components; used for graphics 
;*************************************************
  div = uv2dvF_Wrap(u,v)                ; u,v ==> divergence; (:,:,:,:)

  ud  = new ( dimsizes(u), typeof(u), "No_FillValue")
  vd  = new ( dimsizes(v), typeof(v), "No_FillValue")
  dv2uvf(div,ud,vd)                     ; divergence ==> divergent components  

  copy_VarCoords(u, ud ) 
  copy_VarCoords(u, vd ) 
  ud@long_name  = "Zonal Divergent Wind"
  ud@units      = u@units
  vd@long_name  = "Meridional Divergent Wind"
  vd@units      = v@units

if (PRINT_RESULT) then
    printVarSummary(ud)                 ; (time,level,lat,lon)
    printMinMax(ud,0)
    print("-----")
    printVarSummary(vd)                 ; (time,level,lat,lon)
    printMinMax(vd,0)
    print("-----")
end if

;*************************************************
; plot results
;*************************************************    

  scl5  = 1e5                                  ; arbitrary: used for nicer plot values
  sclab5= "(10~S~-5~N~)"                       ; used later   
  SCLAB5= "(10~S~5~N~)"           

  scl6  = 1e6  
  sclab6= "(10~S~-6~N~)"         
  SCLAB6= "(10~S~6~N~)"         

  plot := new(2,graphic)

  wks   = gsn_open_wks("png","mfc_div")        ; send graphics to PNG file
  resd                 = True
  resd@cnFillOn        = True                  ; color
  resd@cnLinesOn       = False                 ; turn off contour lines

  resd@cnLevelSelectionMode = "ManualLevels"   ; set manual contour levels
  resd@cnMinLevelValF  = -15.                  ; set min contour level
  resd@cnMaxLevelValF  =  15.                  ; set max contour level
  resd@cnLevelSpacingF =   1.                  ; set contour spacing
 ;resd@cnFillPalette   = "cmocean_balance"     ; NCL 6.5.0
  resd@cnFillPalette   = "ViBlGrWhYeOrRe"

  resd@mpFillOn        = False                 ; turn off map fill
  resd@vcRefMagnitudeF = 3.                    ; make vectors larger
  resd@vcRefLengthF    = 0.025                 ; reference vector length
  resd@vcGlyphStyle    = "CurlyVector"         ; turn on curly vectors
  resd@vcMinDistanceF  = 0.010                 ; thin the vectors
  resd@vcRefAnnoOrthogonalPosF = -1.0          ; move ref vector up
  resd@gsnLeftString   = "Divergent Wind"
  resd@gsnScalarContour= True                  ; vectors over contours
  
  LEVP = 700
  DIV  = div(nt,{LEVP},:,:)                     ; keep meta data
  DIV  = DIV*scl6                               ; nicer numbers                 

  resd@tiMainString    = "Divergence and Divergent Winds" 
  resd@gsnCenterString = LEVP+"hPa: "+date     
  resd@gsnRightString  = sclab6+" "+div@units
  dplt = gsn_csm_vector_scalar_map(wks,ud(nt,{LEVP},:,:),vd(nt,{LEVP},:,:),DIV,resd)

;--- Moisture Transport [uq, vq] 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@lbLabelBarOn      = False            ; turn off individual cb's
  res@mpFillOn          = False            ; turn off map fill
                                           ; Use a common scale
  res@cnLevelSelectionMode = "ManualLevels"; manual set levels so lb consistent
  res@cnMaxLevelValF       =  140.0        ; max level
  res@cnMinLevelValF       = -res@cnMaxLevelValF     ; min level
  res@cnLevelSpacingF      =   10.0        ; contour interval

  LEVP    = 700
  res@gsnCenterString      = LEVP+"hPa"
  plot(0) = gsn_csm_contour_map(wks,uq(nt,{LEVP},:,:),res)
  plot(1) = gsn_csm_contour_map(wks,vq(nt,{LEVP},:,:),res)

  resP                     = True                ; modify the panel plot
  resP@gsnPanelMainString  = date+": Unweighted Moisture Flux Components"
  resP@gsnPanelLabelBar    = True                ; add common colorbar
  gsn_panel(wks,plot,(/2,1/),resP)               ; now draw as one plot

;--- Integrated Moisture Transport [iuq, ivq]

  delete(res@gsnCenterString)              ; not used for this plot
  res@cnMaxLevelValF       = 10.0          ; min level
  res@cnMinLevelValF       = -res@cnMaxLevelValF     ; min level
  res@cnLevelSpacingF      =  0.5          ; contour interval

  IUQ     = iuq(nt,:,:)                    ; local array: keep meta data
  IUQ     = IUQ/scl5                       ; scale for plot
  res@gsnRightString  = SCLAB5+" "+iuq@units
  plot(0) = gsn_csm_contour_map(wks,IUQ,res)

  IVQ     = ivq(nt,:,:)                    ; local array: keep meta data
  IVQ     = IVQ/scl5
  res@gsnRightString  = SCLAB5+" "+ivq@units
  plot(1) = gsn_csm_contour_map(wks,IVQ,res)

  resP@gsnPanelMainString  = date+": Mass Wgt. Component Moisture Flux"
  gsn_panel(wks,plot,(/2,1/),resP)               ; now draw as one plot

  delete( [/IUQ, IVQ/] )                   ; no longer needed

;---Divergence of Moisture Flux

  res@cnMaxLevelValF       = 100.0          ; min level
  res@cnMinLevelValF       = -res@cnMaxLevelValF     ; min level
  res@cnLevelSpacingF      =  5.0          ; contour interval

  LEVP    = 700
  DUVQ    = duvq(nt,{LEVP},:,:)                    ; keep meta data
  DUVQ    = DUVQ*scl6                              ; scale for plot
  res@gsnCenterString = LEVP+"hPa"
  res@gsnRightString  = sclab6+" "+duvq@units
  plot(0) = gsn_csm_contour_map(wks,DUVQ,res)

  LEVP    = 500
  DUVQ    = duvq(nt,{LEVP},:,:)                    ; keep meta data
  DUVQ    = DUVQ*scl6
  res@gsnCenterString = LEVP+"hPa"
  res@gsnRightString  = sclab6+" "+duvq@units
  plot(1) = gsn_csm_contour_map(wks,DUVQ,res)

  resP@gsnPanelMainString  = date+": Divergence of Moisture Flux"
  gsn_panel(wks,plot,(/2,1/),resP)                ; now draw as one plot

  delete(DUVQ)                                        ; no longer needed
  delete([/res@gsnCenterString, res@gsnRightString/]) ; not used in next plot 

;---Integrated Divergence of Moisture Flux Convergence [no scaling]

  res@gsnDraw              = True
  res@gsnFrame             = True
  res@lbLabelBarOn         = True        

 ;res@cnFillPalette        = "cmp_flux"
  res@cnMaxLevelValF       =  0.50                ; min level
  res@cnMinLevelValF       = -res@cnMaxLevelValF  ; min level
  res@cnLevelSpacingF      =  0.050               ; contour interval
  res@tiMainString         = date+": VIMFC"

  plt = gsn_csm_contour_map(wks,VIMFC(nt,:,:) ,res)

mfc_div_2.ncl

mfc_div_2.ncl: The above MFC equation can be partitioned as follows:

代码语言:javascript
复制
MFC => Moisture Flux Convergence

MFC_advect = -(u*(dq/dx)+v*(dq/dy))   ; advect moisture term
MFC_conv   =  -q*((du/dx)+(dv/dy) )   ; [con/div]ergence term
MFC        = MFC_advect + MFC_conv

The MFC_advect can be derived using advect_variable for global rectilinear grids or advect_variable_cfd for regional rectilinear grids The MFC_conv can be derived using: uv2dvF_Wrap or uv2dvG_Wrap for global rectilinear grids or uv2dv_cfd for regional rectilinear grids. Then, multiply the derived quantity by specific humidity [q].

代码语言:javascript
复制
;*************************************************
; mfc_div_2.ncl
;
; Similar to mfc_div_1.ncl except a different approach is used.
;
; Rather tha using DIV(U*Q) directly, this script expands this into two separate components.
;
; Concepts illustrated:
;   MFC = Moisture Flux Convergence
;
;   MFC_advect = -(u*(dq/dx)+v*(dq/dy) )    ; advection term
;   MFC_conv   = -q*((du/dx)+(dv/dy) )      ; con(div)-vergence
;
;   MFC = MFC_advect + MFC_convection         
;
;   - Plot a number of quantities
;*************************************************
;---Calculate the Horizontal Moisture Flux Convergence [MFC]
;*************************************************
;---High frequency source data: hourly/3hr/6hr/12hr/daily .... NOT monthly values
;---References:
;---http://www.cgd.ucar.edu/cas/catalog/newbudgets/
;---http://tornado.sfsu.edu/geosciences/classes/e260/AtmosphericRivers/Moisture%20Flux.pdf
;---https://www.spc.noaa.gov/publications/banacos/mfc-sls.pdf
;===================================================================
;   Data Source: ESRL Physical Sciences Division
;        https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html
;   NCEP Reanalysis data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, 
;   from their Web site at https://www.esrl.noaa.gov/psd/
;===================================================================

  ptop = 300              ; 'shum' upper level                  
  ptop@units = "hPa"
  g    = 9.80665          ; m/s2

  date  = 20080715        ; NH summer

;---ESRL: CDC data

  diri = "./"
  filq = "shum.2008.nc"   ; daily data for current year [366 days]
  filu = "uwnd.2008.nc"
  filv = "vwnd.2008.nc"
  filps= "pres.sfc.2008.nc"

  pthu = diri+filu      
  pthv = diri+filv
  pthq = diri+filq     
  pthps= diri+filps   

  fu   = addfile(pthu ,"r")
  fv   = addfile(pthv ,"r")
  fq   = addfile(pthq ,"r")
  fps  = addfile(pthps,"r")

;---Time

  ymd  = cd_calendar(fu->time, -2)    ; ymd[*]: human readable
  nt   = ind(ymd.eq.date)             ; date for plotting and testing
 
TEST = True        
if (.not.TEST) then                   ; all times                   
  u    = fu->uwnd(:,{1000:ptop},:,:)  ; m/s, (time,level,lat,lon)
  v    = fv->vwnd(:,{1000:ptop},:,:)
  q    = fq->shum                     ; [kg/kg], 1000-300 levels only
  ps   = fps->pres                    ; Pa=>[kg/(m-s2)], (time,lat,lon)

else                                  ; one time step; keep time dimension [ nt:nt: ]                     
  u    = fu->uwnd(nt:nt,{1000:ptop},:,:); m/s, (time,level,lat,lon)
  v    = fv->vwnd(nt:nt,{1000:ptop},:,:)
  q    = fq->shum(nt:nt,:,:,:)        ; [kg/kg], 1000-300 levels only
  ps   = fps->pres(nt:nt,:,:)         ; Pa=>[kg/(m-s2)], (time,lat,lon)
  nt   = 0                            ; only one time step
end if

;---Vertical levels

  ptop = ptop*100
  ptop@units = "Pa"

  plev = q&level                      ; hPa  
  plev = plev*100                     ; [100000,...,30000] Pa [kg/(m-s2)]
  plev@units = "Pa"

;---Change [kg/kg] to [g/kg]; not necessary: but common units for q

  q    = q*1000            
  q@units = "g/kg"

;---Divergence function [used later] requires S->N grid order

  u    = u(:,:,::-1,:)  
  v    = v(:,:,::-1,:)
  q    = q(:,:,::-1,:)     
  ps   =ps(:,  ::-1,:)       

;---Layer thickness: ; Pa=>[kg/(m-s2)], (time,level,lat,lon) 
;---Mass weighting:  (dp/g) => [Pa/(m/s2)] => (Pa-s2)/m => [kg/(m-s2)][s2/m] =>  (kg/m2)
;---Reference: http://www.cgd.ucar.edu/cas/catalog/newbudgets/

  dp   = dpres_plevel_Wrap(plev, ps, ptop, 0) ; Pa; layar thickness 

  dpg  = dp/g    
  dpg@long_name = "Layer Mass Weighting"
  dpg@units     = "kg/m2"      ; dp/g, Pa/(m s-2), reduce to kg m-2

;************************************************
; Calculate the MFC_advection term
;   MFC_advect = -(u*(dq/dx)+v*(dq/dy) ) 
; Internally, gradients are calculated via spherical harmonics
;*************************************************    


  long_name = "MFC_advection"
  units     = "g/(kg-s)"       ; (m/s)*(g/kg)*(1/m) => (m/s)*(g/kg-m) => g/(kg-s)

  gridType  = 1   ; global fixed grid ordered S->N
  opt_adv   = 0   ; return only the advected variable; no gradients
  mfc_adv   = advect_variable(u,v,q, gridType, long_name, units, opt_adv)
  mfc_adv   = -mfc_adv

  printVarSummary(mfc_adv)  
  printMinMax(mfc_adv, 0)  
  print("--------")

;************************************************
; Calculate the MFC_convergence term
;   MFC_conv   = -q*((du/dx)+(dv/dy) )      ; con(div)-vergence
;*************************************************    

  duv  = uv2dvF_Wrap(u, v)        ; (1/m)(m/s) => (1/s) ; (time,level,lat,lon)

  mfc_con   = -q*duv           
  mfc_con@long_name = "MFC_convergence"
  mfc_con@units     = "g/(kg-s)"  ; (g/kg)(1/s) => g/(kg-s)
  copy_VarCoords(duv,mfc_con)
  delete(duv)

;************************************************
; Calculate the total MFC
;*************************************************    

  mfc = mfc_adv + mfc_con
  mfc@long_name = "Moisture Flux Convergence"
  mfc@units     = "g/(kg-s)"  ; (g/kg)(1/s) => g/(kg-s)
  

PRINT_RAW = True
if (PRINT_RAW) then
  printVarSummary(mfc_adv)                          ; (time,level,lat,lon)
  printMinMax(mfc_adv,0)
  printVarSummary(mfc_con)                          
  printMinMax(mfc_con,0)
  print("-----")
  printVarSummary(mfc)                          
  printMinMax(mfc,0)
  print("-----")

  printVarSummary(ps)                         ; (time,lat,lon); Pa => kg/(m-s2) 
  printMinMax(ps,0)
  print("-----")
  printVarSummary(dp)                         ; (time,level,lat,lon); Pa => kg/(m-s2)
  printMinMax(dp,0)
  print("-----")
                                ; examine layer thickness at selected locations
  print(dp(nt,:,{40},{180}))    ; mid-Pacific
  print(dp(nt,:,{40},{255}))    ; Boulder, CO 
  print("-----")
end if

;---Integrated mass weighted moisture flux components

  mfc_adv_dpg = mfc_adv*dpg                ; mass weighted 'uq'; [m/s][g/kg][kg/m2]=>[m/s][g/kg]
  imfc_adv    =  dim_sum_n(mfc_adv_dpg, 1)
  imfc_adv@long_name = "Integrated Mass Flux Advection" 
  imfc_adv@LONG_NAME = "Sum: Mass Weighted Integrated Mass Flux Advection: mfc_adv*dpg" 
  imfc_adv@units     = "[m/s][g/kg]"
  copy_VarCoords(u(:,0,:,:), imfc_adv); (time,lat,lon)
  delete(mfc_adv_dpg)

  mfc_con_dpg = mfc_con*dpg                ; mass weighted 'mfc_con'; [m/s][g/kg][kg/m2]=>[m/s][g/kg] 
  imfc_con    = dim_sum_n(mfc_con_dpg, 1)
  imfc_con@long_name = "Integrated  Mass Flux Convergence"
  imfc_con@LONG_NAME = "Sum: Mass Weighted Integrated Mass Flux Convergence [mfc_con*dpg]" 
  imfc_con@units     = "[m/s][g/kg]"
  copy_VarCoords(v(:,0,:,:), imfc_con); (time,lat,lon)
  delete(mfc_con_dpg)

  VIMFC =  imfc_adv +  imfc_con    
  VIMFC@long_name = "VIMFC"
  VIMFC@LONG_NAME = "VIMFC: [imfc_adv+imfc_con]"
  copy_VarCoords(q(:,0,:,:),VIMFC)            ; (time,lat,lon)

PRINT_RESULT = True
if (PRINT_RESULT) then
    printVarSummary(imfc_adv)                 ; (time,lat,lon)
    printMinMax(imfc_adv,0)
    print("-----")
    printVarSummary(imfc_con)                 ; (time,lat,lon)
    printMinMax(imfc_con,0)
    print("-----")
    printVarSummary(VIMFC)               ; (time,lat,lon)
    printMinMax(VIMFC,0)
    print("-----")
end if

;*************************************************
; plot results
;*************************************************    

  scl5  = 1e5                                  ; arbitrary: used for nicer plot values
  sclab5= "(10~S~-5~N~)"                       ; used later   
  SCLAB5= "(10~S~5~N~)"           

  scl6  = 1e6  
  sclab6= "(10~S~-6~N~)"         
  SCLAB6= "(10~S~6~N~)"         

  plot := new(2,graphic)

  wks   = gsn_open_wks("png","mfc_div_2")        ; 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@mpFillOn          = False            ; turn off map fill
  res@lbLabelBarOn      = False            ; turn off individual cb's
                                           ; Use a common scale
  res@cnLevelSelectionMode = "ManualLevels"; manual set levels so lb consistent
  res@cnMaxLevelValF       =   12.0        ; max level
  res@cnMinLevelValF       = -res@cnMaxLevelValF     ; min level
  res@cnLevelSpacingF      =    0.5        ; contour interval

  LEVP    = 700
  res@gsnCenterString      = LEVP+"hPa"

  MFC_ADV = mfc_adv(nt,{LEVP},:,:)         ; keep meta data
  MFC_ADV = MFC_ADV*scl5
  res@gsnRightString  = sclab5+" "+mfc_adv@units
  plot(0) = gsn_csm_contour_map(wks,MFC_ADV,res)

  MFC_CON = mfc_con(nt,{LEVP},:,:)
  MFC_CON = MFC_CON*scl5
  res@gsnRightString  = sclab5+" "+mfc_con@units
  plot(1) = gsn_csm_contour_map(wks,MFC_CON,res)

  resP                     = True                ; modify the panel plot
  resP@gsnPanelMainString  = date+": Unweighted MFC_ADV, MFC_CON"
  resP@gsnPanelLabelBar    = True                ; add common colorbar
  gsn_panel(wks,plot,(/2,1/),resP)               ; now draw as one plot

;--- Integrated Moisture Transport [iuq, ivq]

  delete([/res@gsnRightString, res@gsnCenterString/]) ; not used for this plot
    
  res@cnMaxLevelValF       =  0.50                   ; min level
  res@cnMinLevelValF       = -res@cnMaxLevelValF     ; min level
  res@cnLevelSpacingF      =  0.05                   ; contour interval

  IMFC_ADV = imfc_adv(nt,:,:)                    ; local array: keep meta data
  plot(0)  = gsn_csm_contour_map(wks,IMFC_ADV,res)

  IMFC_CON = imfc_con(nt,:,:)                    ; local array: keep meta data
  plot(1) = gsn_csm_contour_map(wks,IMFC_CON,res)

  resP@gsnPanelMainString  = date+": Integrated Moisture Flux: Advect, Convergence"
  gsn_panel(wks,plot,(/2,1/),resP)               ; now draw as one plot

  delete( [/IMFC_ADV, IMFC_CON/] )                   ; no longer needed

  res@lbLabelBarOn      = True   
  res@gsnDraw = True
  res@gsnFrame= True

;---Integrated Divergence of Moisture Flux Convergence [no scaling]

 ;res@cnFillPalette        = "cmp_flux"
  res@cnLevelSelectionMode = "ManualLevels"; manual set levels so lb consistent
  res@cnMaxLevelValF       =  0.50                ; min level
  res@cnMinLevelValF       = -res@cnMaxLevelValF  ; min level
  res@cnLevelSpacingF      =  0.050               ; contour interval
  res@tiMainString         = date+": VIMFC: [IMFC_ADV+IMFC_CON]"

  plt = gsn_csm_contour_map(wks,VIMFC(nt,:,:) ,res)

本文参与 腾讯云自媒体分享计划,分享自微信公众号。
原始发表:2020-12-17,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 MeteoAI 微信公众号,前往查看

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

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

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • mfc_div_1.ncl
  • mfc_div_2.ncl
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档