;*******************************************************************************************
;
; Script to interpolate data at hybrid model levels into chosen constant pressure levels:
;
;   - Uses "vinth2p_ecmwf" to interpolate into pressure levels, with extrapolation
;     into pressure leves below the surface pressure. An ECMWF formulation is used
;     for the extrapolation of temperature and geopotential data.
;
;*******************************************************************************************

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"

begin

;;; »»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»» DATA FILES TO INTERPOLATE ««««««««««««««««««««««««««««««««««««««««««
; ============================================= TEMPERATURE =====================================================
Ta = addfile("ta_6hrLev_MPI-ESM1-2-LR_historical_r1i1p1f1_gn_197901.nc", "r")

; =================== VERTICAL WIND ======================================================================
;Wa = addfile("wap_6hrLev_MPI-ESM1-2-LR_historical_r1i1p1f1_gn_197001.nc", "r")

; =================== GEOPOTENTIAL HEIGHT ================================================================
Zg = addfile("zg_6hrLev_MPI-ESM1-2-LR_historical_r1i1p1f1_gn_197901.nc", "r")

; ========================================= SURFACE GEOPOTENTIAL ================================================
Zs = addfile("orog_fx_MPI-ESM1-2-LR_historical_r1i1p1f1_gn.nc", "r")





; --------------------------------------- Constants ----------------------------------------------------
p0 = 1000.0     ; Surface reference pressure (hPa).

; The normal acceleration of gravity:
; g  = 9.80616 - This value has been chosen as the most representative of the acceleration of gravity
;                at sea-level at latitude 45º in accordance with the opinion of the International
;                Association of Geodesy (1950). This is the value generally used by physicists.
;
; g  = 9.80665 - Value adopted by WMO (1966) which agrees with that recommended by the International
;                Committee of Weights and Measures (1901, 1948).

g  = 9.80665

; ------------------------------------------------------------------------------------------------------

; ################# LOADING THE NEEDED VARIABLES ###############################################
phis    = Zs->orog*g          ; surface geopotential (m^2/sec^2)


time = Ta->time      ; time (time:calendar = "proleptic_gregorian")

lat  = Ta->lat           ; latitudes

lon  = Ta->lon         ; longitudes

ps   = Ta->ps           ; surface pressures in Pa.

hyAM = Ta->ap       ; hybrid A coefficients at layer midpoints ordered bottom-to-top.
                                ; These are expected to be normalized to one (1.0). If not, divide by p0 (Pa).

hyBM = Ta->b         ; hybrid B coefficients at layer midpoints ordered bottom-to-top.

; Hybrid coefficients ordered top-to-bottom.
hyam = hyAM(::-1)/(p0*100)
hybm = hyBM(::-1)




;; The 23 pressure levels to interpolate (the same as ERA-40):
;lev = (/1,2,3,5,7,10,20,30,50,70,100,150,200,250,300,400,500,600,700,775,850,925,1000/)
lev = (/1,5,10,20,30,50,70,100,150,200,250,300,400,500,600,700,850,925,1000/);;; 19-plevs as in the "Amom" table


lev!0                  = "lev"             ; variable/dim name 
lev&lev              =  lev               ; create coordinate variable
lev@long_name = "pressure"   ; attach some attributes
lev@units          = "hPa"
lev@positive      = "down"

                             ; 1 = linear
intyp = 2              ; 2 = log (Variables assumed to vary with logarithm of pressure)
                            ; 3 = log-log

; Get dimension sizes:
ntim = dimsizes(time)
nlev = dimsizes(lev)
nlat = dimsizes(lat)
nlon = dimsizes(lon)

; »»»»»»»»»»»»»»»»»»»»»»» TEMPERATURE ««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««

;;; Getting temperature at lowest model level

Th     = Ta->ta(:,::-1,:,:)   ; air temperature (K) at hybrid levels
                                        ; ordered top-to-bottom and south-to-north.

; Assign named dimensions
Th!0 = "time"
Th!1 = "lev"
Th!2 = "lat"
Th!3 = "lon"

Dm     = dimsizes(Th)             ; Get number of dimensions (time,lev,lat,lon)

tb     = Th(:,Dm(1)-1,:,:)        ; temperature at the lowest level.

delete(Th)





; »»»»»»»»»»»»»»»»»»»»»»» GEOPOTENTIAL HEIGHT ««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««
varflg = -1                 ; geopotential height is variable

kxtrp  = True             ; Logical variable:
                                  ; False = no extrapolation when the pressure level is outside of
                                  ; the range of surface pressure. True = extrapolate.

Zh     = Zg->zg(:,::-1,:,:)   ; geopotential height at hybrid levels reversed top-to-bottom.


; Geopotential height interpolated into pressure levels with extrapolation below the surface pressure:
zg = vinth2p_ecmwf(Zh,hyam,hybm,lev,ps,intyp,p0,1,kxtrp,varflg,tb,phis)

delete(Zh)


;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
; Writting output to netCDF file with associated meta data:
;++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

;------ The output file name in native grid:
        Zout = "zg_6hrPLev_MPI-ESM1-2-LR_historical_r1i1p1f1_gn_197901.nc"

	system("/bin/rm -f " + Zout)   ; remove if exists

	fZout  = addfile (Zout, "c")   ; open output file

;==================================================================================================
; explicitly declare file definition mode. Improve efficiency.
;==================================================================================================
        setfileoption(fZout,"DefineMode",True)

;==================================================================================================
; Copy file attributes:
;==================================================================================================

;------ Assign file attributes:
    fAtt                                         = True
	fAtt@institution                     = "Max Planck Institute for Meteorology"
	fAtt@institute_id                    = "MPI-M"
    fAtt@experiment_id               = "historical"
    fAtt@parent_variant_label     = "r1i1p1f1"
	fAtt@parent_experiment_id   = "piControl"
    fAtt@experiment                    = "historical"
	fAtt@frequency                      = "6hr"
    fAtt@project_id                      = "CMIP6"
	fAtt@table_id                         = "Table 6hrLev"
    fAtt@modeling_realm           = "atmos"
	fAtt@realization                    = "1"

	fileattdef( fZout, fAtt )

;==================================================================================================
; predefine the coordinate variables and their dimensionality.
; To get an UNLIMITED record dimension, the dimensionality is set to -1 and the unlimited array to True.
;==================================================================================================
        dimNames = (/"time", "lev", "lat", "lon"/)
	dimSizes = (/ -1   , nlev ,  nlat, nlon /)
	dimUnlim = (/ True , False, False, False/)
	filedimdef(fZout,dimNames,dimSizes,dimUnlim)

;==================================================================================================
; Defining the variable's dimension name(s). 
;==================================================================================================
       filevardef(fZout, "time"          ,typeof(time)    ,"time")
       filevardef(fZout, "lev"           ,typeof(lev)        ,"lev")
       filevardef(fZout, "lat"           ,typeof(lat)         ,"lat")
       filevardef(fZout, "lon"           ,typeof(lon)       ,"lon")
       filevardef(fZout, "zg"           ,typeof(zg)          ,(/"time","lev","lat","lon"/))

;==================================================================================================
; Copy attributes associated with each variable to the file
; All attributes associated with each variable will be copied.
;==================================================================================================
       filevarattdef(fZout,"zg",zg)            ; copy zg attributes
       filevarattdef(fZout,"time",time)          ; copy time attributes
       filevarattdef(fZout,"lev",lev)            ; copy lev attributes
       filevarattdef(fZout,"lat",lat)            ; copy lat attributes
       filevarattdef(fZout,"lon",lon)            ; copy lon attributes

;===================================================================
; explicitly exit file definition mode. **NOT REQUIRED**
;===================================================================
        setfileoption(fZout,"DefineMode",False)

;===================================================================
; output only the data values since the dimensionality and such have
; been predefined. The "(/", "/)" syntax tells NCL to only output the
; data values to the predefined locations on the file.
;====================================================================

       fZout->time = (/time/)
       fZout->lev  = (/lev/)
       fZout->lat  = (/lat/)
       fZout->lon  = (/lon/)
       fZout->zg  = (/zg/)

; »»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»»X««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««««





; Deleting variables with a time dimension, because it will change in the next iteration (i.e. month):
delete(time)
delete(zg)
delete(tb)
delete(ps)
delete(lev)


end