#!/usr/bin/magml
# Metview Macro

#  **************************** LICENSE START ***********************************
# 
#  Copyright 2013 ECMWF. This software is distributed under the terms
#  of the Apache License version 2.0. In applying this license, ECMWF does not
#  waive the privileges and immunities granted to it by virtue of its status as
#  an Intergovernmental Organization or submit itself to any jurisdiction.
# 
#  ***************************** LICENSE END ************************************
# 
#---------------------------------------
# Define global parameters
#---------------------------------------

global VDF_FILE=""
global OUTPUT_DIR=""
global WORK_DIR=""
global LOG_FILE=""

global STEP_NUMBER=-1
global SURF_PARAMS=""
global UPPER_PARAMS=""
global VERTICAL_GRID=""
global ELEVATION_PARAM=""
global MIN_VERTICAL="0"
global MAX_VERTICAL="0"
global VDF_LEVEL=""

global AREA=""
global GRID=""

global INPUT_GRB=nil
global DATE_LST=nil
global LEV_LST=nil
global LEV_TYPE=""
global LON_WEST=0
global LON_EAST=0
global LAT_NORTH=0
global LAT_SOUTH=0
global NX=0
global NY=0
global NZ=0
global REVERT_LAT=0
global REVERT_LON=0
global errH=0
global infoH=0

#--------------------------------
# The vapor setup command
#--------------------------------

global VAPOR_SETUP_CMD=". $MV_VAPOR_BIN/vapor-setup.sh"

#----------------------------------------------------
#The exe to convert grib into raw binary for vapor
#----------------------------------------------------

global TORAW_EXE=getenv("METVIEW_BIN") & "/gribToVapor"

#---------------------------------------------------------------
# Read arguments
#---------------------------------------------------------------

#Get arguments
global ARGS
ARGS=arguments()
global CURRENT_ARG_INDEX=1

if count(ARGS) = 0 then
	log_fail("No arguments specified!")
end if

print("Arguments: ",ARGS)

#Get output dir 
OUTPUT_DIR=ARGS[1]
WORK_DIR=OUTPUT_DIR

#Get vdf filename
VDF_FILE=nextArg() & ".vdf"

#Get vdf check mode
val=nextArg()
doCheck=0
if val = "1" then
    doCheck=1
    if exist(VDF_FILE) then
    	log_info("VDF file already exists! No conversion is performed!")
    	stop()
    end if	
end if 

#Get vdf refinement level
VDF_LEVEL=nextArg()

#Get vertical grid params
VERTICAL_GRID=nextArg()
if VERTICAL_GRID = "LAYERED" then
     ELEVATION_PARAM=nextArg()
	 MIN_VERTICAL=nextArg()
	 MAX_VERTICAL=nextArg()
end if

#Get interpoltaion params
val=nextArg()
doInterpolate=0
if val = "1" then
    doInterpolate=1
    AREA=parse(nextArg(),"/")
    GRID=parse(nextArg(),"/")
end if 

#Get step number
STEP_NUMBER=number(nextArg())

#Get parameter list
SURF_PARAMS=parse(nextArg(),"/")
UPPER_PARAMS=parse(nextArg(),"/") 

#Get working dir
WORK_DIR=nextArg()

#Get input files
inputMode=nextArg()
if inputMode = "ICON" then
    INPUT_GRB=nil
	while hasNextArg() do 
			INPUT_GRB=INPUT_GRB & read(nextArg())
	end while		
else
	log_fail("InputMode=" & inputMode & " not implemented!")
end if	

#---------------------------------------
#Log input parameters
#---------------------------------------

log_info("Input arguments --->")
log_info(tab & "VDF_FILE=" & VDF_FILE)
log_info(tab & "OUTPUT_DIR=" & OUTPUT_DIR)
log_info(tab & "SURF_PARAMS=" &  listToString(SURF_PARAMS))
log_info(tab & "UPPER_PARAMS=" & listToString(SURF_PARAMS))
log_info(tab & "VERTICAL_GRID=" & VERTICAL_GRID)
if VERTICAL_GRID = "LAYERED" then
     log_info(tab & "ELEVATION_PARAM=" & ELEVATION_PARAM)
	 log_info(tab & "MIN_VERTICAL=" & MIN_VERTICAL)
	 log_info(tab & "MAX_VERTICAL=" & MAX_VERTICAL)
end if
if doInterpolate = 1 then
    log_info(tab & "Interpolate to AREA=" & listToString(AREA) & " GRID=" & listToString(GRID))
end if
log_info(tab & "STEP_NUMBER=" & STEP_NUMBER)
log_info(tab & "VDF_LEVEL=" & VDF_LEVEL)

#---------------------------------------
# Check directories
#---------------------------------------

if not(exist(OUTPUT_DIR)) then
   if shell("mkdir -p " & OUTPUT_DIR) <> 0 then
      print("Could not create output directory!")
      print("Shell command failed: mkdir -p " & OUTPUT_DIR)
      fail()
   end if   
end if

if not(exist(WORK_DIR)) then
   if shell("mkdir -p " & WORK_DIR) <> 0 then
      print("Could not create working directory!")
      print("Shell command failed: mkdir -p " & WORK_DIR)
      fail()
   end if   
end if

#------------------------------------
# Interpolation
#------------------------------------

#Should be done only for the required fields
if doInterpolate = 1 then	 
  	INPUT_GRB=read(
         	data : INPUT_GRB, 
		 	area: AREA,
		 	grid: GRID       
       	)
end if
  
#------------------------------------
# Get vertical structure 
#------------------------------------
getVertical()

#------------------------------------
# Get dates and times
#------------------------------------
getDateTime()

#------------------------------------
# Get geography
#------------------------------------
getGeography()

#------------------------------------
# Sorting
#------------------------------------

if LEV_TYPE = "ml" or LEV_TYPE = "pl" then
   INPUT_GRB=sort(INPUT_GRB,["param","date","time","step","levelist"],["<","<","<","<",">"])
#   #Ig=sort(INPUT_GRB,"parameter")
#   #INPUT_GRB=sort(INPUT_GRB)
else
   INPUT_GRB=sort(INPUT_GRB,["param","date","time","step","levelist"],["<","<","<","<","<"])
end if

#------------------------------------
# Check vertical parameters
#------------------------------------

if VERTICAL_GRID = "LAYERED" then
	prepareElevation()
end if	

#-------------------------
# Create the vdf file
#-------------------------
createVdf()  

#-------------------------
# Process 2D params
#-------------------------

log_info("Processing 2D params --->")

loop v in SURF_PARAMS 
  
  log_info(tab & "param: " & v)
  
  g=read(
         	data : INPUT_GRB, 
		 	levtype: "sfc",
		 	param: v		       
       	)
  if g =nil then
  	   log_fail("No data found for param: " & v) 	   
  end if	
  
  varName=checkVarName(v,"2D") 
  processVar(g,varName,"2D")
      
end loop   

#-------------------------
#Process 3D params
#-------------------------
log_info("Processing 3D params --->")

loop v in UPPER_PARAMS 
  
  log_info(tab & "param: " & v)
  
  g=read(
         data : INPUT_GRB, 
		 levtype: LEV_TYPE,
		 param: v		       
       )
  if g =nil then
  	   log_fail("No data found for param: " & v) 	   
  end if	
        
  varName=checkVarName(v,"3D")      
  processVar(g,varName,"3D")
    
end loop   

stop(0)

#======== Functions ===================================

#=======================================================
# Read the all the validity dates for the first available
# variable. The result is stored in DATE_LST
#=======================================================

function getDateTime()

 log_info("Check dates --->")
 
 par=SURF_PARAMS[1]
 g=read(data : INPUT_GRB, param: par)
 
 DATE_LST=nil
 loop fg in g	
        
     if STEP_NUMBER = -1 or  count(DATE_LST) <  STEP_NUMBER then
    
     	vd=grib_get_long(fg,"validityDate")
     	vt=grib_get_long(fg,"validityTime")
    
     	d=(date: vd,time: vt)        
        
     	if count(DATE_LST) =0 then
     		DATE_LST = DATE_LST & [d]
 	 	else 
 	    	found = 0
 	 		loop item in DATE_LST
 	 	   		if item["date"] = d["date"] and 
 	 	      		item["time"] = d["time"] then	 	     
 	 	          		found = 1 
 	 	   		end if 
 	 		end loop 
 	 	
 	 		if found = 0 then
 	 			DATE_LST = DATE_LST & [d]
 	 		end if        
 	 	end if
 	 end if
 	 		
 end loop		
 
 log_info(tab &"Validity dates:")
 for i=1 to count(DATE_LST) do
 	log_info(tab & i & ") date=" & DATE_LST[i].date & ", time=" & DATE_LST[i].time)
 end for	
 
 #print("dates")
 #print(DATE_LST)
 
end getDateTime

#=======================================================
# Read the all the level related parameters from the first
# upper level parameter.
#=======================================================

function getVertical()

  log_info("Check vertical --->")
  
  par=UPPER_PARAMS[1]
  if par = "lnsp" then
     if count(UPPER_PARAMS) < 2 then
  		log_fail("No upper level params specified!")
  	 end if	
     par=UPPER_PARAMS[2]
  end if
     
  g=read(data : INPUT_GRB, param: par)
  levType="sfc"
  cnt=1
  while levType = "sfc" and cnt <> count(g) do
  	  levType=grib_get_string(g[cnt],"levelType")
  	  cnt=cnt+1
  end while
  
  if levType = "sfc" then
  	 log_fail("No vertical level is found for param: " & par)
  end if	
  
  #Set the level type
  LEV_TYPE=levType
 
  #Generate the level list 
  loop fg in g	
     lev=grib_get_string(fg,"level")
     if count(LEV_LST) = 0 then
     	LEV_LST=LEV_LST & [lev]
     else if lev not in LEV_LST then
        LEV_LST=LEV_LST & [lev] 	
     end if
  end loop
  
  NZ=count(LEV_LST)
   
  log_info(tab & "Level type: " & LEV_TYPE)
  log_info(tab & "Level num: " & NZ)
  levStr=""
  for i=1 to count(LEV_LST) do
  	if i = 1 then
  		levStr=levStr & LEV_LST[i]
  	else
  		levStr=levStr & "," & LEV_LST[i]
  	end if
  end for
  log_info(tab & "Levels: " & levStr)
  
  			
  
end getVertical

function getGeography()

	log_info("Check geography ---> ")
	 
	g=INPUT_GRB[1]
	
	lat1=grib_get_double(g,"latitudeOfFirstGridPointInDegrees")
	lat2=grib_get_double(g,"latitudeOfLastGridPointInDegrees")
	
	if lat1 > lat2 then
		LAT_NORTH=lat1
		LAT_SOUTH=lat2
	else
		LAT_NORTH=lat2
		LAT_SOUTH=lat1	
	end if
		
	lon1=grib_get_double(g,"longitudeOfFirstGridPointInDegrees")
	lon2=grib_get_double(g,"longitudeOfLastGridPointInDegrees")
	
	#log_info(tab & "lon1/lon2: " & lon1 & "/" &lon2)
		
	if lon1 > 180 then
		lon1 = lon1-360
	end if	
	if lon2 > 180 then
		lon2 = lon2-360
	end if
	
	#log_info(tab & "lon1/lon2: " & lon1 & "/" &lon2)

	if lon1 < lon2 then
		LON_WEST=lon1
		LON_EAST=lon2
	else
		LON_WEST=lon2
		LON_EAST=lon1
	end if
		
	NX=grib_get_long(g,"Ni")
	NY=grib_get_long(g,"Nj")
	
	log_info(tab & "nx/ny: " & NX & "/" & NY)
	log_info(tab & "N/W/S/E: " & LAT_NORTH & "/" & LON_WEST & "/" & LAT_SOUTH & "/" & LON_EAST)	
    #log_info(tab & "revert lat: " & REVERT_LAT)
    #log_info(tab & "revert lon: " & REVERT_LON)

end getGeography
	
#====================================================
# Check and modify an original parameter names to 
# fit it to VAPOR's needs
#====================================================		

function checkVarName(varName,varType:string)

  if varName = "z" then  
    if varType = "3D" then
  		return "ELEVATION"
	else if varType = "2D" then
		return "HGT"
	end if
  end if
  	
  numberLst=["0","1","2","3","4","5","6","7","8","9"]
  if type(varName) = "number" then
  		return "v" & string(varName)
  else		
  	s = substring(varName,1,1)
  	if s in numberLst then
     	return "v" & varName
    end if
  end if  
  
  return varName
end checkVarName     

#====================================================
# Create the VDF file
#====================================================

function createVdf() 
	
  #Write al the timestamps into a ascii file.
  #vdfcreate will need it!	  
  fName=WORK_DIR & "/tsfile.txt"
  f=file(fName) 
  for i=1 to count(DATE_LST) do  
    vdate=DATE_LST[i]["date"]
    vtime=DATE_LST[i]["time"]	   
    write(f,vdate &  vtime, newline)
  end for    
  f=0  
  
  #Generate 2d parameter string for vdfcreate
  par2D=""
  for i=1 to count(SURF_PARAMS) do
	if i=1 then
		par2D=checkVarName(SURF_PARAMS[i],"2D")
	else
		par2D=par2D & ":" & checkVarName(SURF_PARAMS[i],"2D")	
	end if
  end for
	
  #Generate 2d parameter string for vdfcreate
  par3D=""
  for i=1 to count(UPPER_PARAMS) do
	if i=1 then
		par3D=checkVarName(UPPER_PARAMS[i],"3D")
	else
		par3D=par3D & ":" & checkVarName(UPPER_PARAMS[i],"3D")	
	end if
  end for
	
  fLog=WORK_DIR & "/log.txt"
  
  if VERTICAL_GRID = "LAYERED" then
     gridType="layered"
  else
  	 gridType="regular"
  	 if LEV_TYPE = "ml" or LEV_TYPE = "pl" then
  	     MIN_VERTICAL="-" & LEV_LST[1]
  	     MAX_VERTICAL="-" & LEV_LST[count(LEV_LST)]
  	 else
  	     MIN_VERTICAL=LEV_LST[1]
  	     MAX_VERTICAL=LEV_LST[count(LEV_LST)]
  	 end if    
  end if 	    
  
  log_info("---> Create VDF file with vdfcreate!")
  
  #Create VDF file
  geoFactor=111177
  cmd=VAPOR_SETUP_CMD & ";" &
      " cd " & OUTPUT_DIR & ";" & 
      " vdfcreate  -dimension " & NX & "x" & NY & "x" & NZ &
          " -gridtype " & gridType &
          " -mapprojection \"+proj=latlong  +ellps=sphere\" " &
          " -level " & VDF_LEVEL &
          " -numts " & count(DATE_LST) &
 #         " -usertimes tsfile.txt " &
          " -extents " & LON_WEST*geoFactor & ":" &
                LAT_SOUTH*geoFactor & ":" &
               	MIN_VERTICAL & ":" &
               	LON_EAST*geoFactor & ":" &
                LAT_NORTH*geoFactor & ":" &
                MAX_VERTICAL & 
          " -vars3d " & par3D &
          " -vars2dxy " & par2D &
          " " & VDF_FILE &  
          " > " & fLog            

  print(cmd)
  if shell(cmd) <> 0 then
  	  log_error("vdfcreate failed!")
  	  log_error("command: " & cmd)
  	  #log_error("output " & read(fLog))
  	  fail()
  else
  	log_info("command: " & cmd)		
  end if
  	
end createVdf

function toGrib(g:fieldset,ts:number) 
	
	#Save grib file
    grbFile=WORK_DIR & "/tmp_" & ts & ".grib"
    write(grbFile,g)

end toGrib    
    

#====================================================
# Write one timestep of one parameter to vdf/vdc
#====================================================

function toVdf(g:fieldset,varName:string,ts:number,varType:string) 
	
	#Save grib file
    #grbFile=WORK_DIR & "/" & varType & "_" & varName & "_" & ts & ".grib"
    #grbFile=WORK_DIR & "/tmp_" & ts & ".grib"
    grbFile=WORK_DIR & "/tmp.grib"
    write(grbFile,g)
    
    #log_info(" grib: " & varName & " " & count(g))
    
    mv_bin=getenv("METVIEW_BIN")
    
    #Generate the raw binary file
    #rawFile=WORK_DIR & "/" & varType & "_" & varName & "_" & ts & ".raw"
    #rawFile=WORK_DIR & "/tmp_" & ts & ".raw"
    rawFile=WORK_DIR & "/tmp.raw"
    cmd=TORAW_EXE & " " & grbFile & " " & rawFile
    
    ret= shell(cmd)
    if ret <> 0 then
  	   log_error("Conversion of GRIB to raw binary failed!")
  	   log_fail("command: " & cmd)
    end if		
   
	#Add raw binary to vdf/vdc
    cmd=VAPOR_SETUP_CMD & ";" &
        " cd " & OUTPUT_DIR & ";" & 
        " raw2vdf -ts " & ts &
                " -varname " & varName &
                " " & VDF_FILE &
                " " & rawFile
 
    ret= shell(cmd)
    if ret <> 0 then
  	   log_error("Adding raw binary to VAPOR failed!")
  	   log_fail("command: " &cmd)
    end if	                 

end toVdf  
     
function processVar(g:fieldset,varName:string,varType:string)     
       
  for i=1 to count(DATE_LST) do
    
    vdate=DATE_LST[i]["date"]
    vtime=DATE_LST[i]["time"]	
    
    log_info(tab & "--> processing: " & vdate & " " & vtime)
     	
    gRes=nil 	
    loop f in g       
      if grib_get_long(f,"validityDate") = vdate and
         grib_get_long(f,"validityTime") = vtime then
         gRes=gRes & f
      end if
    end loop        
    
    #Write to vdf
    toVdf(gRes,varName,i-1,varType)
    
   end for 
   
end processVar   
 
function prepareElevation()

	log_info("Checking vertical parameters --->")
	
	g=read(data : INPUT_GRB, 
	       param: ELEVATION_PARAM,
	       levtype: LEV_TYPE)
	
	if g = nil then
		#log_fail("Elevation parameter=" &  ELEVATION_PARAM & "not found in input data!")
		if LEV_TYPE  = "ml" and ELEVATION_PARAM = "z" then
		    log_info(tab & "Geopotential is not present! Need to be computed or model levels")
		    computeElevation_ml()
		else
			log_fail("Elevation parameter=" &  ELEVATION_PARAM & " not found in input data!") 
		end if	   
	end if	      
	    
	if ELEVATION_PARAM not in UPPER_PARAMS then
       UPPER_PARAMS = UPPER_PARAMS & [ELEVATION_PARAM]
    end if   
	       
end prepareElevation 
 
function computeElevation_ml()
	
	 t=read(data : INPUT_GRB, 
	       param: "t",
	       levtype: LEV_TYPE)
	    
	 q=read(data : INPUT_GRB, 
	       param: "q",
	       levtype: LEV_TYPE)
	    
	 lnsp=read(data : INPUT_GRB, 
	       param: "lnsp",
	       levtype: LEV_TYPE,
	       level: 1)
	       
	 zs=read(data : INPUT_GRB, 
	       param: "z",
	       levtype: "sfc")      
	       	 
	 if t=nil then 
	 	log_fail("Parameter=t not found on model levels! Height of model levels cannot be computed!")
	 else if q = nil then 
	    log_fail("Parameter=q not found on model levels! Height of model levels cannot be computed!")	 	 
	 else if lnsp = nil then 
	    log_fail("Parameter=lnsp not found! Height of model levels cannot be computed!") 	 
	 else if zs = nil then 
	    log_fail("Surface geopotential is not found! Height of model levels cannot be computed!")
	 end if
	 	 	 
	 gRes=nil 
	 for i=1 to count(DATE_LST) do
	 	
	 	vdate=DATE_LST[i]["date"]
    	vtime=DATE_LST[i]["time"]
 
        tAct=nil 	
        loop f in t        
             if grib_get_long(f,"validityDate") = vdate and
                grib_get_long(f,"validityTime") = vtime then
                tAct=tAct & f
             end if
        end loop
        
        qAct=nil 	
        loop f in q        
             if grib_get_long(f,"validityDate") = vdate and
                grib_get_long(f,"validityTime") = vtime then
                qAct=qAct & f
             end if
        end loop      
          
        lnspAct=nil 	
        loop f in lnsp        
             if grib_get_long(f,"validityDate") = vdate and
                grib_get_long(f,"validityTime") = vtime then
                lnspAct=lnspAct & f
             end if
        end loop 
        
        zsAct=nil 	
        loop f in zs        
             if grib_get_long(f,"validityDate") = vdate and
                grib_get_long(f,"validityTime") = vtime then
                zsAct=zsAct & f
             end if
        end loop 
                       
 	    gRes=gRes & geopotential_on_ml(tAct,qAct,lnspAct,zsAct)
 
    end for

    INPUT_GRB = INPUT_GRB & gRes
    
end computeElevation_ml

 
function nextArg()
  
  CURRENT_ARG_INDEX=CURRENT_ARG_INDEX+1   
  
  if CURRENT_ARG_INDEX <= count(ARGS) then
  	    return ARGS[CURRENT_ARG_INDEX]
  else
  	    return ""
  end if
  	    
end nextArg  

function hasNextArg()
	if CURRENT_ARG_INDEX < count(ARGS) then 
	  return 1
	else 
      return 0
    end if
end hasNextArg         
	
function log_fail(txt:string)
	
	log_error(txt)
	fail()
	
end log_fail	
	
function log_error(txt:string)
	
	if errH=0 then 
		errH=file(WORK_DIR & "/m_err.txt")
	end if		
	write(errH,txt,newline)

end log_error	

function log_info(txt:string)

	if infoH=0 then 
		infoH=file(WORK_DIR & "/m_log.txt")
	end if
	write(infoH,txt,newline)

end log_info


function geopotential_on_ml (T_fs:fieldset, q_fs:fieldset, lnsp_fs:fieldset, zs_fs:fieldset)

    Rd = 287.06
    g  = 9.80665

    #Do some basic checking on the input data - 
    #only check the first field of each fielset for efficiency
    gridtype_t = grib_get_string(T_fs[1],    'gridType')
    gridtype_q = grib_get_string(q_fs[1],    'gridType')
    gridtype_l = grib_get_string(lnsp_fs[1], 'gridType')
    gridtype_z = grib_get_string(zs_fs[1],   'gridType')

    if (gridtype_t = 'sh' or gridtype_q = 'sh' or gridtype_l = 'sh' or gridtype_z = 'sh') then
        log_fail("geopotential_on_ml: fields must be gridded, not spectral")
    end if

    #How many levels are we computing?
    levelSize = count(T_fs) 
	bottomLevel = grib_get_long(T_fs[1], 'level')
    topLevel = grib_get_long(T_fs[levelSize], 'level')

    if bottomLevel < topLevel then
       log_fail("geopotential_on_ml: fields are not sorted in descending order according to model levels")
    end if 
    
    #print("levels: ",bottomLevel," ",topLevel)

    zs   = values(zs_fs)     # extract the data for this field as a vector
    lnsp = values(lnsp_fs)   # extract the data for this field as a vector

    pv = grib_get_double_array(lnsp_fs, 'pv') # for computing pressures at model levels
    sp = exp(lnsp)                          # surface pressure
    z_h = zs                                # orography

    #Get the coefficients for computing the pressures  
    #e.g. L-91: A=1..92, B=93..184
    pv = grib_get_double_array(lnsp_fs, 'pv')
    totalLevNum=count(pv)/2-1
    A = pv[1, totalLevNum+1]
    B = pv[totalLevNum+2, 2*totalLevNum+2]  
    
    if totalLevNum <> bottomLevel then
       log_fail("geopotential_on_ml: bottom level should be " & totalLeveNum & " instead of " & bottomLevel)
    end if   
      
    #A = pv[1, levelSize+1]                #e.g. L-91: A=1..92, B=93..184
    #B = pv[levelSize+2, levelSize+levelSize+2]

    z_f = z_h        # output vector
    z_out_fs = nil   # output fieldset - the end result

    #Compute the bottom pressure (on half-levels)
    #initialise to the lowest numbered level
    Ph_levplusone = (A[totalLevNum+1] + (B[totalLevNum+1] * sp)) 

    # We want to integrate up into the atmosphere, starting at the ground
    # so we start at the lowest level (highest number) and keep
    # accumulating the height as we go.
    # See the IFS documentation:
    # http://www.ecmwf.int/research/ifsdocs/DYNAMICS/Chap2_Discretization4.html
    # For speed and file I/O, we perform the computations with vectors instead
    # of fieldsets.

    for lev=bottomLevel to topLevel by -1 do

		levIdx=bottomLevel-lev+1 #levelSize-lev+1
		
		print("level: ",lev," ",levIdx)
		
        # we assume that the data are in (numerically) descending level order - check if true

        T_level = grib_get_long(T_fs[levIdx], 'level')
        q_level = grib_get_long(q_fs[levIdx], 'level')

        if (T_level <> lev) then
            log_fail ('T field index ', levIdx, ' should be level ', lev, ' but it is ', T_level)
        else if (q_level <> lev) then
            log_fail ('q field index ', levIdx, ' should be level ', lev, ' but it is ', q_level)
        end if


        q = values(q_fs[levIdx]) # extract the data for this field as a vector
        T = values(T_fs[levIdx]) # extract the data for this field as a vector

        # quick check for data consistency

        if ((count(q) <> count(T)) or
            (count(q) <> count(sp)) or
            (count(q) <> count(z_h))) then
            errmsg = 'mvl_geopotential_on_ml: T, Q, LNSP and Z, must have the same number of grid points (they have ' &
                      count(T) & ', ' & count(q) & ', ' & count(sp) & ' and ' & count(z_h) & ' respectively).'
            log_fail(errmsg)
        end if

        T = T*(1.+0.609133*q)   # compute moist temperature
        q = 0                   # free the memory for 'q'

        # compute the pressures (on half-levels)

        Ph_lev = (A[lev] + (B[lev] * sp))

        if lev = 1 then
            dlogP = log(Ph_levplusone/0.1) # Ph[lev] is zero, so don't divide by it
            alpha = log(2)
        else
            dlogP = log(Ph_levplusone/Ph_lev)
            dP    = Ph_levplusone-Ph_lev
            alpha = 1 - Ph_lev/dP*dlogP
        end if

        TRd = T*Rd
        T   = 0      # free memory for 'T'

        # z_f is the geopotential of this full level
        # integrate from previous (lower) half-level z_h to the full level

        z_f = z_h + (TRd*alpha)

        # z_h is the geopotential of 'half-levels'
        # integrate z_h to next half level

        z_h = z_h + (TRd*dlogP)

        # store the result (z_f) in a field and add to the output fieldset
        # (add it to the front, not the end, because we are going 'backwards'
        # through the fields)
        
        out_field = set_values(T_fs[levIdx], z_f) 
        out_field = grib_set_string(out_field, ['shortName', "z"])       
        out_field = grib_set_long(out_field,  ['generatingProcessIdentifier',128,
                                              'level', lev])
             
        z_out_fs = z_out_fs &out_field

        Ph_levplusone  = Ph_lev    # store for the next iteration, will be equivalent to (A[lev+1] + (B[lev+1] * sp))

    end for

    return z_out_fs

end geopotential_on_ml

function listToString(lst)

  s=""
  for i=1 to count(lst) do
  	if i = 1 then
  		s=s & lst[i]
  	else
  		s=s & "," & lst[i]
  	end if
  end for
  
  return s
  
end listToString