source: configs/dcmip2016/PLOTS/testcase/BAROCLINIC_WAVE/dcmip1-lat-lon-terminator_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
;
; Modified by Kevin Viner for terminator ouput
; 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"   
; ================================================;

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 - average integrated atomic Cl
  ;                     2 - average integrated diatomic Cl
  ;                     3 - average integrated Cly
        do var_choice = 1,3

        if (var_choice .eq. 1) then                             ; Select Cl
                filename_string = "DYNAMICO-case161latlonQ1int"         ; output using X11
        else if (var_choice .eq. 2) then                        ; Select Cl2
                filename_string = "DYNAMICO-case161latlonQ2int"         ; output using png
        else if (var_choice .eq. 3) then                        ; Select Cly
                filename_string = "DYNAMICO-case161latlonQtot_int"
        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

        output_per_day=8 ; 3-hourly outputs
        do d = 1,30*output_per_day, (output_per_day/4)
        day_curr = d/output_per_day

        varload1 = f->Q1_col_int
        varload2 = f->Q2_col_int
        if (var_choice .eq. 1) then                             ; Select Cl
                var1 = varload1(d,:,:)/4e-6
                filename_string = "DYNAMICO-case161latlonQ1int"         ; output using X11
        else if (var_choice .eq. 2) then                        ; Select Cl2
                var1 = varload2(d,:,:)/4e-6
                filename_string = "DYNAMICO-case161latlonQ2int"         ; output using png
        else if (var_choice .eq. 3) then                                ; Select Cly
                var1 = (varload1(d,:,:)+2*varload2(d,:,:))/4e-6 ; should equal 1.
                print("Max Cly error :")
                print(max(abs(var1-1.)))
                var1 = 1e5*abs(var1-1.)
                filename_string = "DYNAMICO-case161latlonQtot_int"
        end if
        end if
        end if


        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, Integrated Average Cl"                     ; Plot title if required
        pres@txString = pltTitle

                ; Change the output type and name

        gsn_define_colormap(wks,"gui_default")  

        res1@cnMaxLevelValF  = 1.2                              ; max contour color label
        res1@cnMinLevelValF  = -0.1                             ; min contour color label
        res1@cnLevelSpacingF = 0.1                              ; contour color spacing                 ; choose a colormap
 
        else if (var_choice .eq. 2) then        ;================ Large-scale precipitation plot ====================
 
        pltTitle="DYNAMICO Case 161, Integrated Average Cl2"    ; 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  = 0.48                     ; max contour color label
        res1@cnMinLevelValF  = -0.04                            ; min contour color label
        res1@cnLevelSpacingF = 0.04                             ; contour color spacing                 ; choose a colormap
 
        else if (var_choice .eq. 3) then        ;================ Surface Pressure plot ====================
 
        pltTitle="DYNAMICO Case 161, Integrated Average Cly"                    ; Plot title if required
        pres@txString = pltTitle


        var1@units = " "
;       var2@units = " "
;       var3@units = " "
;       var4@units = " "
                ; Change the output type and name

        gsn_define_colormap(wks,"gui_default")  

        res1@cnMaxLevelValF  = 1.24                             ; max contour color label
        res1@cnMinLevelValF  = 0.78                             ; min contour color label
        res1@cnLevelSpacingF = 0.04                             ; 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.52,txres3)

        txres3@txAngleF       = 0.
        gsn_text_ndc(wks,"Longitude",0.54,0.28,txres3)

;       draw(plot)
        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