source: configs/dcmip2016/PLOTS/testcase/BAROCLINIC_WAVE/dcmip1-lat-lon-TPLSPS_LOOP.ncl @ 437

;================================================;
;       Example ncl script to produce the set of
;       lat-lon plots for DCMIP-2016 test case 1
; Created by James Kent, Christiane Jablonowski
;       and Paul Ullrich (University of Michigan) for DCMIP-2012 
;
; Modified by Colin Zarzycki for DCMIP-2016
; Further modified by Maximo Menchaca for creating gifs
;
; v1.01 - 6/7/16 
; v1.02 - 6/8/16 - Typo in var_choice options
;================================================;
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"   
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"   
; ================================================;

; NOTE, if your model does not output T850 OR P at level
; midpoints OR you cannot calculate P from hybrid levels, you
; will need to make model specific additions denoted by
; CALCPHERE below

; PRECL needs to be in m/s for contours to plot correctly

begin

;=================================================;
; open file and read in data 
; GIVEN AS (time,lev,lat,lon) from 0 to n-1
;=================================================;

  ; NOTE, that this file needs to be a catted file including all times at daily interval
  ; if you are outputting 1 file per time, you can run "ncrcat dcmip1.files.*.nc cat.nc"
  ; for one file

        f = addfile("output_dcmip2016_regular.nc","r")

use_pdf = False

                ; Input useful parameters

        lat  = f->lat
        lon  = f->lon
        lev  = f->lev
        nlat   = getfilevardimsizes(f, "lat" )
        nlon   = getfilevardimsizes(f, "lon" )
        nlev   = getfilevardimsizes(f, "lev" )

  ; We want plots of temperature, vertical velocity and relative vorticity 
  ; at the 850 hPa level. We also want surface pressure.

  ; Select var_choice   1 - Temperature at 850 hPa
  ;                     2 - Large-scale precipitation
  ;                     3 - Surface Pressure

        do var_choice = 1,3
;       var_choice = 1 
        
        if (var_choice .eq. 1) then                             ; Select T

                if (isfilevar(f,"T850")) then

                        varload = f->T850

                else    

                        ; we interp T to the 850 hPa level
                        varload0 = f->T                         ; T is a 4D var
                        varload  = f->T(:,0,:,:)
                        ; Calculate P using hybrid coefficients
                        ; If you have pressure as a variable just do
                        ; P = f->P
      ; CALCPHERE
                        hyam = f->hyam
                hybm = f->hybm
                P0   = f->P0
                        PS   = f->PS
                        P    = 0.0*varload0
                        P = pres_hybrid_ccm(PS,P0,hyam,hybm)
                        P!0 = "time"
                        P!1 = "lev"
                        P!2 = "lat"
                        P!3 = "lon"
                        ; interp to 850 hPa level
                        plevel = 85000.0
                        varload = int2p_n(P,varload0,plevel,2,1)

                end if

        filename_string = "DYNAMICO-case161latlonT850"

        else if (var_choice .eq. 2) then                        ; Select Large-scale precipitation

                varload = f->PRECT
    varload=varload*8.64e7              
        filename_string = "DYNAMICO-case161latlonPRECT"

        else if (var_choice .eq. 3) then                        ; Select PS
        filename_string = "DYNAMICO-case161latlonPS"

                varload = f->PS

        end if
        end if
        end if

if(use_pdf) then
        wks  = gsn_open_wks("pdf",filename_string)      ; output using X11
else
        wks  = gsn_open_wks("png",filename_string)      ; output using png
end if

                ; We want the output at days 6, 9, 12 and 15. If your model 
                ; output is in 6 hourly intervals this becomes indices 24, 36, 48 and 60

        do d = 1,240,2
          day_curr = d/8

;       day6  = 24 
;       day9  = 36
;       day12 = 48
;       day15 = 60

        var1 = varload(d,:,:)
;       var2 = varload(day9,:,:)
;       var3 = varload(day12,:,:)
;       var4 = varload(day15,:,:)

                ; We don't want long-name in the plots

        var1@long_name = " "
;       var2@long_name = " "
;       var3@long_name = " "
;       var4@long_name = " "

                ; Delete loaded data

;       delete(varload)

                ; We now produce the plot

;       plot  = new (4, graphic)                                ; define plot - need 4 panels

        res1                      = True
        res1@gsnDraw              = False                       ; panel plot 
        res1@gsnFrame             = False                       ; don't draw yet
        res1@cnFillOn             = True
        res1@cnLinesOn            = True
        res1@gsnSpreadColors      = True
        res1@lbLabelAutoStride    = True
        res1@gsnCenterString      = ""
        res1@tiMainString         = ""
        res1@vpWidthF             = 0.38
        res1@vpHeightF            = 0.19
        res1@cnLevelSelectionMode = "ManualLevels"
        res1@cnInfoLabelOn        = False                       ; don't give each
        res1@cnLineLabelsOn       = False                       ; panel its own
        res1@lbLabelBarOn         = False                       ; label bar

        pres                          = True
        pres@gsnMaximize              = True  
        pres@gsnPanelLabelBar         = True                    ; Communal label bar
        pres@gsnPanelLeft             = 0.1
        pres@gsnPanelRight            = 0.9
        pres@pmLabelBarOrthogonalPosF = -0.03
        pres@gsnFrame                 = False
        pres@lbLabelStride            = 1

        res1@sfXArray        = lon                              ; uses lon as plot x-axis
        res1@sfYArray        = lat                              ; uses lat for y axis
        res1@trYReverse      = False                            ; reverses y-axis, false
        res1@tiYAxisString   = ""                               ; y-axis title
        res1@tiXAxisString   = ""                               ; x-axis title

        if (var_choice .eq. 1) then     ;================ Temperature plot ====================
 
        pltTitle="DYNAMICO Case 161, T850 hPa"                  ; Plot title if required
        pres@txString = pltTitle

                ; Change the output type and name

        gsn_define_colormap(wks,"gui_default")  

        res1@cnMaxLevelValF  = 300.0                            ; max contour color label
        res1@cnMinLevelValF  = 230.0                            ; min contour color label
        res1@cnLevelSpacingF = 10.0                             ; contour color spacing                 ; choose a colormap
 
        else if (var_choice .eq. 2) then        ;================ Large-scale precipitation plot ====================
 
        pltTitle="DYNAMICO Case 161, Large-scale precipitation"         ; Plot title if required
        pres@txString = pltTitle

        var1@units = " "
;       var2@units = " "
;       var3@units = " "
;       var4@units = " "

        res1@cnLinesOn            = False                               ; Turn lines off for clearer plot

                ; Change the output type and name

        gsn_define_colormap(wks,"gui_default")  

                ; Note that the maximum might be larger than these contour spacing values

        res1@cnMaxLevelValF  = 50.0                             ; max contour color label
        res1@cnMinLevelValF  = 0.0                              ; min contour color label
        res1@cnLevelSpacingF = 4.0                              ; contour color spacing                 ; choose a colormap
 
        else if (var_choice .eq. 3) then        ;================ Surface Pressure plot ====================
 
        pltTitle="DYNAMICO Case 161, PS"                        ; Plot title if required
        pres@txString = pltTitle

                ; Convert to hPa 

        var1=var1/100.0
;       var2=var2/100.0
;       var3=var3/100.0
;       var4=var4/100.0

        var1@units = "hPa "
;       var2@units = "hPa "
;       var3@units = "hPa "
;       var4@units = "hPa "

                ; Change the output type and name

        gsn_define_colormap(wks,"gui_default")  

        res1@cnMaxLevelValF  = 1020.0                           ; max contour color label
        res1@cnMinLevelValF  = 920.0                            ; min contour color label
        res1@cnLevelSpacingF = 10.0                             ; contour color spacing                 ; choose a colormap

        end if                          ;============== end plot choice ===============
        end if                          ;============== end plot choice ===============
        end if                          ;============== end plot choice ===============

        res1@gsnCenterString = "t = " + day_curr + " days"
        plot = gsn_csm_contour(wks,var1(:,:),res1)              ; plot var1

;       res1@gsnCenterString = "t = 9"
;       plot(1) = gsn_csm_contour(wks,var2(:,:),res1)           ; plot var2
;
;       res1@gsnCenterString = "t = 12"
;       plot(2) = gsn_csm_contour(wks,var3(:,:),res1)           ; plot var3
;
;       res1@gsnCenterString = "t = 15"
;       plot(3) = gsn_csm_contour(wks,var4(:,:),res1)           ; plot var4

        gsn_panel(wks,plot,(/1,1/),pres)                        ; 2x2 plot

                ; Add latitude and longitude labels

        txres3                = True
        txres3@txAngleF       = 90.
        txres3@txFontHeightF  = 0.02
        gsn_text_ndc(wks,"Latitude",0.08,0.49,txres3)

        txres3@txAngleF       = 0.
        gsn_text_ndc(wks,"Longitude",0.5,0.22,txres3)

        frame(wks)

        end do
        delete(varload)
 cmd = "convert -delay 25 "+ filename_string+"*.png "+filename_string+".gif"
 system(cmd)
  cmd = "rm "+filename_string+"*.png"
 system(cmd)
        end do

end