Index: flexpart_code/fpmetbinary_mod.F90
===================================================================
--- flexpart_code/fpmetbinary_mod.F90	(revision fd86dea2732396822196778365c05ed2ef15f881)
+++ flexpart_code/fpmetbinary_mod.F90	(revision 8624a75e20c0cc6dff42b4eb97182a138f9aa1e7)
@@ -238,16 +238,43 @@
             ncret = nf90_def_dim(ncid, 'nxmax', nxmax, nxmax_dimid)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","maximum dimension of wind fields in x")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
             ncret = nf90_def_dim(ncid, 'nymax', nymax, nymax_dimid)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","maximum dimension of wind fields in y")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
             ncret = nf90_def_dim(ncid, 'nzmax', nzmax, nzmax_dimid)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","maximum dimension of wind fields in z")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
             ncret = nf90_def_dim(ncid, 'nuvzmax', nuvzmax, nuvzmax_dimid)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description"," maximum dimension of (u,v) wind field in z")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
             ncret = nf90_def_dim(ncid, 'nwzmax', nwzmax, nwzmax_dimid)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description"," maximum dimension of w component of  wind field in z")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
             ncret = nf90_def_dim(ncid, 'maxspec', maxspec, maxspec_dimid)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description"," maximum number of species")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
             ncret = nf90_def_dim(ncid, 'numclass', numclass, numclass_dimid)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description"," maximum number of ageclasses")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ! There are a handful of variables indexed from 0 to n, rather than 0 to n-1,
@@ -255,19 +282,36 @@
             ncret = nf90_def_dim(ncid, 'zero_to_nzmax', nzmax+1, zero_to_nzmax_dimid)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","variable to change indexing of variables with nzmax dimension")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
             ncret = nf90_def_dim(ncid, 'zero_to_maxnests', maxnests+1, zero_to_maxnests_dimid)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","variable to change indexing of variables with maxnests dimension")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ! This is for a couple of small arrays that store polar stereographic stuff
             ncret = nf90_def_dim(ncid, 'polemap_dim', 9, polemap_dimid)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","variable to store polar stereographic indexing")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ! These two values come from conv_mod
             ncret = nf90_def_dim(ncid, 'nconvlevmax_dim', nconvlevmax, nconvlevmax_dimid)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","maximum number of levels for convection")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
             ncret = nf90_def_dim(ncid, 'na_dim', na, na_dimid)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","maximum number of levels for convection +1")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ! Scalar values
-
 
             dim1dids = (/preproc_fmt_str_dimid/)
@@ -288,4 +332,7 @@
             ncret = nf90_put_var(ncid, ncvarid, nx)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in x")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'ny', NF90_INT, ncvarid)
@@ -293,4 +340,7 @@
             ncret = nf90_put_var(ncid, ncvarid, ny)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in y")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'nxmin1', NF90_INT, ncvarid)
@@ -298,4 +348,7 @@
             ncret = nf90_put_var(ncid, ncvarid, nxmin1)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in x -1 (nx-1)")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'nymin1', NF90_INT, ncvarid)
@@ -303,4 +356,8 @@
             ncret = nf90_put_var(ncid, ncvarid, nymin1)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in y -1 (ny-1)")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'nxfield', NF90_INT, ncvarid)
@@ -308,4 +365,7 @@
             ncret = nf90_put_var(ncid, ncvarid, nxfield)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in x for limited area fields")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'nuvz', NF90_INT, ncvarid)
@@ -313,4 +373,7 @@
             ncret = nf90_put_var(ncid, ncvarid, nuvz)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of (u,v) wind fields in z direction")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'nwz', NF90_INT, ncvarid)
@@ -318,4 +381,7 @@
             ncret = nf90_put_var(ncid, ncvarid, nwz)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of w wind fields in z direction")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'nz', NF90_INT, ncvarid)
@@ -323,4 +389,7 @@
             ncret = nf90_put_var(ncid, ncvarid, nz)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","number of vertical levels in the transformed coordinates")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'nmixz', NF90_INT, ncvarid)
@@ -328,4 +397,7 @@
             ncret = nf90_put_var(ncid, ncvarid, nmixz)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","number of levels up to maximum PBL height (3500 m)")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'nlev_ec', NF90_INT, ncvarid)
@@ -333,4 +405,7 @@
             ncret = nf90_put_var(ncid, ncvarid, nlev_ec)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","number of vertical levels ecmwf model")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'dx', NF90_FLOAT, ncvarid)
@@ -338,4 +413,7 @@
             ncret = nf90_put_var(ncid, ncvarid, dx)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","grid distance in x direction")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'dy', NF90_FLOAT, ncvarid)
@@ -343,4 +421,7 @@
             ncret = nf90_put_var(ncid, ncvarid, dy)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","grid distance in y direction")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'xlon0', NF90_FLOAT, ncvarid)
@@ -348,4 +429,7 @@
             ncret = nf90_put_var(ncid, ncvarid, xlon0)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","geographical longitude of the lower left corner of the wind fields")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'ylat0', NF90_FLOAT, ncvarid)
@@ -353,4 +437,7 @@
             ncret = nf90_put_var(ncid, ncvarid, ylat0)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","geographical latitude of the lower left corner of the wind fields")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'dxconst', NF90_FLOAT, ncvarid)
@@ -358,4 +445,7 @@
             ncret = nf90_put_var(ncid, ncvarid, dxconst)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable needed for utransform")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'dyconst', NF90_FLOAT, ncvarid)
@@ -363,6 +453,7 @@
             ncret = nf90_put_var(ncid, ncvarid, dyconst)
             call handle_nf90_err(ncret)
-
-
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable needed for vtransform")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ! Fixed fields, static in time
@@ -380,4 +471,8 @@
 &                                oro(0:nxmax-1, 0:nymax-1))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","orography of the ECMWF model")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
 
             ncret = nf90_def_var(ncid, 'excessoro', NF90_FLOAT, &
@@ -392,4 +487,8 @@
 &                                excessoro(0:nxmax-1, 0:nymax-1))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","excess orography mother domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'lsm', NF90_FLOAT, &
@@ -404,4 +503,8 @@
 &                                lsm(0:nxmax-1, 0:nymax-1))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","land sea mask")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             dim3dids = (/nxmax_dimid, nymax_dimid, numclass_dimid/)
@@ -418,4 +521,8 @@
 &                                xlanduse(0:nxmax-1, 0:nymax-1, 1:numclass))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","area fraction in percent")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             dim1dids = (/nzmax_dimid/)
@@ -431,5 +538,7 @@
 &                                height(1:nzmax))
             call handle_nf90_err(ncret)
-
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","heights of all model levels")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m")
 
 
@@ -449,4 +558,7 @@
 &                                uu(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","U component of wind in the [horizontal] direction")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
 
             ncret = nf90_def_var(ncid, 'vv', NF90_FLOAT, &
@@ -461,4 +573,8 @@
 &                                vv(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","V component of wind in the Y[horizontal] direction")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
 
             ncret = nf90_def_var(ncid, 'uupol', NF90_FLOAT, &
@@ -473,4 +589,8 @@
 &                                uupol(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","horizontal component in polar stereographic projection of wind")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
 
             ncret = nf90_def_var(ncid, 'vvpol', NF90_FLOAT, &
@@ -485,4 +605,8 @@
 &                                vvpol(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","horizontal component in polar stereographic projection of wind")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
 
             ncret = nf90_def_var(ncid, 'ww', NF90_FLOAT, &
@@ -497,4 +621,8 @@
 &                                ww(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","wind component in the Z[vertical] direction")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
 
             ncret = nf90_def_var(ncid, 'tt', NF90_FLOAT, &
@@ -509,4 +637,8 @@
 &                                tt(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","temperature")
+            ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
 
             ncret = nf90_def_var(ncid, 'qv', NF90_FLOAT, &
@@ -521,4 +653,7 @@
 &                                qv(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","specific humidity data")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'pv', NF90_FLOAT, &
@@ -533,4 +668,8 @@
 &                                pv(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","potential vorticity")
+            ncret = nf90_put_att(ncid, ncvarid, "units","K*m**2 kg**-1 s**-1")
+
 
             ncret = nf90_def_var(ncid, 'rho', NF90_FLOAT, &
@@ -545,4 +684,8 @@
 &                                rho(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","air density")
+            ncret = nf90_put_att(ncid, ncvarid, "units","kg m**-3")
+
 
             ncret = nf90_def_var(ncid, 'drhodz', NF90_FLOAT, &
@@ -557,4 +700,8 @@
 &                                drhodz(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","vertical air density gradient")
+            ncret = nf90_put_att(ncid, ncvarid, "units","kg m**-2")
+
 
             ncret = nf90_def_var(ncid, 'clouds', NF90_BYTE, &
@@ -569,5 +716,10 @@
 &                                clouds(0:nxmax-1, 0:nymax-1, 1:nzmax, cm_index))
             call handle_nf90_err(ncret)
-
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","vcloud mask, &
+                                 no cloud no precip = 0, cloud no precip = 1, &
+                                 rainout  conv/lsp dominated  2/3, &
+                                 washout  conv/lsp dominated = 4/5 ")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
 
@@ -586,4 +738,8 @@
 &                                tth(0:nxmax-1, 0:nymax-1, 1:nuvzmax, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","temperature in the original input model level")
+            ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
 
             ncret = nf90_def_var(ncid, 'qvh', NF90_FLOAT, &
@@ -598,4 +754,8 @@
 &                                qvh(0:nxmax-1, 0:nymax-1, 1:nuvzmax, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","specific humidity in the original input model level")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'pplev', NF90_FLOAT, &
@@ -610,4 +770,7 @@
 &                                pplev(0:nxmax-1, 0:nymax-1, 1:nuvzmax, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","pressure in pressure levels for GFS version")
+            ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
 
 
@@ -624,4 +787,7 @@
 &                                cloudsh(0:nxmax-1, 0:nymax-1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","height of the cloud layer")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m")
 
 
@@ -649,4 +815,8 @@
 &                                ps(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","surface pressure")
+            ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
 
             ncret = nf90_def_var(ncid, 'sd', NF90_FLOAT, &
@@ -661,4 +831,8 @@
 &                                sd(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","snow depth")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
 
             ncret = nf90_def_var(ncid, 'msl', NF90_FLOAT, &
@@ -673,4 +847,8 @@
 &                                msl(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","mean sea level pressure")
+            ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
 
             ncret = nf90_def_var(ncid, 'tcc', NF90_FLOAT, &
@@ -685,4 +863,8 @@
 &                                tcc(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","total cloud cover")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'u10', NF90_FLOAT, &
@@ -697,4 +879,7 @@
 &                                u10(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","10 m u component of wind velocity")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
 
             ncret = nf90_def_var(ncid, 'v10', NF90_FLOAT, &
@@ -709,4 +894,8 @@
 &                                v10(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","10 m v component of wind velocity")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
 
             ncret = nf90_def_var(ncid, 'tt2', NF90_FLOAT, &
@@ -721,4 +910,8 @@
 &                                tt2(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","2 m temperature")
+            ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
 
             ncret = nf90_def_var(ncid, 'td2', NF90_FLOAT, &
@@ -733,4 +926,8 @@
 &                                td2(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","2 m dew point temperature")
+            ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
 
             ncret = nf90_def_var(ncid, 'lsprec', NF90_FLOAT, &
@@ -745,4 +942,8 @@
 &                                lsprec(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","large scale total precipitation")
+            ncret = nf90_put_att(ncid, ncvarid, "units","mm h**-1")
+
 
             ncret = nf90_def_var(ncid, 'convprec', NF90_FLOAT, &
@@ -757,4 +958,8 @@
 &                                convprec(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","convective precipitation")
+            ncret = nf90_put_att(ncid, ncvarid, "units","mm h**-1")
+
 
             ncret = nf90_def_var(ncid, 'sshf', NF90_FLOAT, &
@@ -769,4 +974,8 @@
 &                                sshf(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","surface sensible heat flux")
+            ncret = nf90_put_att(ncid, ncvarid, "units","J m**-2")
+
 
             ncret = nf90_def_var(ncid, 'ssr', NF90_FLOAT, &
@@ -781,4 +990,8 @@
 &                                ssr(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","surface solar radiation")
+            ncret = nf90_put_att(ncid, ncvarid, "units","J m**-2")
+
 
             ncret = nf90_def_var(ncid, 'surfstr', NF90_FLOAT, &
@@ -793,4 +1006,8 @@
 &                                surfstr(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","surface stress")
+            ncret = nf90_put_att(ncid, ncvarid, "units","N m**-2 s")
+
 
             ncret = nf90_def_var(ncid, 'ustar', NF90_FLOAT, &
@@ -805,4 +1022,8 @@
 &                                ustar(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","friction velocity")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
 
             ncret = nf90_def_var(ncid, 'wstar', NF90_FLOAT, &
@@ -817,4 +1038,8 @@
 &                                wstar(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","convective velocity scale")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
 
             ncret = nf90_def_var(ncid, 'hmix', NF90_FLOAT, &
@@ -829,4 +1054,8 @@
 &                                hmix(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","mixing height")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
 
             ncret = nf90_def_var(ncid, 'tropopause', NF90_FLOAT, &
@@ -841,4 +1070,8 @@
 &                                tropopause(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","altitude of thermal tropopause")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
 
             ncret = nf90_def_var(ncid, 'oli', NF90_FLOAT, &
@@ -853,4 +1086,8 @@
 &                                oli(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","inverse Obukhov length (1/L)")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
 
             ncret = nf90_def_var(ncid, 'diffk', NF90_FLOAT, &
@@ -865,5 +1102,7 @@
 &                                diffk(0:nxmax-1, 0:nymax-1, 1, cm_index))
             call handle_nf90_err(ncret)
-
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","diffusion coefficient at reference height")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m**2 s**-1")
 
 
@@ -888,5 +1127,7 @@
 &                                vdep(0:nxmax-1, 0:nymax-1, 1:maxspec, cm_index))
             call handle_nf90_err(ncret)
-
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","deposition velocity")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-2")
 
 
@@ -902,4 +1143,7 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                z0(1:numclass))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","surface roughness length")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m")
 
 
@@ -917,4 +1161,8 @@
 &                                akm(1:nwzmax))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","ECMWF vertical discretization parameter")
+            ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
 
             ncret = nf90_def_var(ncid, 'bkm', NF90_FLOAT, &
@@ -929,4 +1177,7 @@
 &                                bkm(1:nwzmax))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","ECMWF  vertical discretization parameter")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
 
@@ -944,4 +1195,8 @@
 &                                akz(1:nuvzmax))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","ECMWF vertical discretization parameter at center of layer")
+            ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
 
             ncret = nf90_def_var(ncid, 'bkz', NF90_FLOAT, &
@@ -956,4 +1211,7 @@
 &                                bkz(1:nuvzmax))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","ECMWF vertical discretization parameter at center of layer")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
 
@@ -971,4 +1229,7 @@
 &                                aknew(1:nzmax))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","model discretization coefficient at the interpolated levels")
+            ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
 
             ncret = nf90_def_var(ncid, 'bknew', NF90_FLOAT, &
@@ -983,4 +1244,7 @@
 &                                bknew(1:nzmax))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","model discretization coefficient at the interpolated levels")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
 
@@ -1017,4 +1281,7 @@
 &                                nxn(1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in x for the nexted domain")
+
 
             ncret = nf90_def_var(ncid, 'nyn', NF90_INT, &
@@ -1029,4 +1296,8 @@
 &                                nyn(1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","actual dimensions of wind fields in y for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'dxn', NF90_FLOAT, &
@@ -1041,4 +1312,8 @@
 &                                dxn(1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","grid distance in x direction for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'dyn', NF90_FLOAT, &
@@ -1053,4 +1328,8 @@
 &                                dyn(1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","grid distance in y direction for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'xlon0n', NF90_FLOAT, &
@@ -1065,4 +1344,9 @@
 &                                xlon0n(1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","geographical longitude of the &
+                                lower left corner of the nested wind fields")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'ylat0n', NF90_FLOAT, &
@@ -1077,6 +1361,8 @@
 &                                ylat0n(1:maxnests))
             call handle_nf90_err(ncret)
-
-
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","geographical latitude of &
+                                 the lower left corner of the nested wind fields")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
 
@@ -1095,4 +1381,8 @@
 &                                oron(0:nxmaxn-1, 0:nymaxn-1, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","orography of the ECMWF model for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
 
             ncret = nf90_def_var(ncid, 'excessoron', NF90_FLOAT, &
@@ -1107,4 +1397,8 @@
 &                                excessoron(0:nxmaxn-1, 0:nymaxn-1, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","excess orography nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'lsmn', NF90_FLOAT, &
@@ -1119,4 +1413,8 @@
 &                                lsmn(0:nxmaxn-1, 0:nymaxn-1, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","land sea mask for the nested domian")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             dim4dids = (/nxmaxn_dimid, nymaxn_dimid, numclass_dimid, maxnests_dimid/)
@@ -1133,4 +1431,8 @@
 &                                xlandusen(0:nxmaxn-1, 0:nymaxn-1, 1:numclass, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","area fraction in percent for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
 !            PRINT *, 'SUM(oron): ', SUM(oron)
@@ -1152,4 +1454,9 @@
 &                                uun(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","U component of wind in the [horizontal] &
+                                 direction for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
 
             ncret = nf90_def_var(ncid, 'vvn', NF90_FLOAT, &
@@ -1164,4 +1471,8 @@
 &                                vvn(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","V component of wind in the Y[horizontal] &
+                                 direction for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
 
             ncret = nf90_def_var(ncid, 'wwn', NF90_FLOAT, &
@@ -1176,4 +1487,9 @@
 &                                wwn(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","wind componend in the Z[vertical] &
+                                  direction for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
 
             ncret = nf90_def_var(ncid, 'ttn', NF90_FLOAT, &
@@ -1188,4 +1504,8 @@
 &                                ttn(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","temperature for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
 
             ncret = nf90_def_var(ncid, 'qvn', NF90_FLOAT, &
@@ -1200,4 +1520,7 @@
 &                                qvn(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","specific humidity data for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'pvn', NF90_FLOAT, &
@@ -1212,4 +1535,7 @@
 &                                pvn(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","potential vorticity for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","K*m**2 kg**-1 s**-1")
 
             ncret = nf90_def_var(ncid, 'rhon', NF90_FLOAT, &
@@ -1224,4 +1550,8 @@
 &                                rhon(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","air density for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","kg m**-3")
+
 
             ncret = nf90_def_var(ncid, 'drhodzn', NF90_FLOAT, &
@@ -1236,4 +1566,7 @@
 &                                drhodzn(0:nxmaxn-1, 0:nymaxn-1, 1:nzmax, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","vertical air density gradient for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","kg m**-2")
 
 
@@ -1252,4 +1585,9 @@
 &                                tthn(0:nxmaxn-1, 0:nymaxn-1, 1:nuvzmax, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","temperature in the original &
+                                 input model level for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
 
             ncret = nf90_def_var(ncid, 'qvhn', NF90_FLOAT, &
@@ -1264,4 +1602,9 @@
 &                                qvhn(0:nxmaxn-1, 0:nymaxn-1, 1:nuvzmax, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","specific humidity in the original &
+                                input model level for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ! Note the new dimensions
@@ -1279,4 +1622,11 @@
 &                                cloudsn(0:nxmaxn-1, 0:nymaxn-1, 0:nzmax, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","vcloud mask, &
+                                 no cloud no precip = 0, cloud no precip = 1, &
+                                 rainout  conv/lsp dominated  2/3, &
+                                 washout  conv/lsp dominated = 4/5 for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ! Note the new dimensions
@@ -1294,4 +1644,8 @@
 &                                cloudsnh(0:nxmaxn-1, 0:nymaxn-1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","cloud mask at the horizontal levels")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
 
@@ -1316,4 +1670,8 @@
 &                                psn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","surface pressure of the nested domian")
+            ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
 
             ncret = nf90_def_var(ncid, 'sdn', NF90_FLOAT, &
@@ -1328,4 +1686,8 @@
 &                                sdn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","snow depth for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m")
+
 
             ncret = nf90_def_var(ncid, 'msln', NF90_FLOAT, &
@@ -1340,4 +1702,8 @@
 &                                msln(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","mean sea level pressure for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
 
             ncret = nf90_def_var(ncid, 'tccn', NF90_FLOAT, &
@@ -1352,4 +1718,8 @@
 &                                tccn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","total cloud cover for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'u10n', NF90_FLOAT, &
@@ -1364,4 +1734,7 @@
 &                                u10n(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","10 m u component of wind velocity for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
 
             ncret = nf90_def_var(ncid, 'v10n', NF90_FLOAT, &
@@ -1376,4 +1749,7 @@
 &                                v10n(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","10 m v component of wind velocity for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
 
             ncret = nf90_def_var(ncid, 'tt2n', NF90_FLOAT, &
@@ -1388,4 +1764,8 @@
 &                                tt2n(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","2 m temperature for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
 
             ncret = nf90_def_var(ncid, 'td2n', NF90_FLOAT, &
@@ -1400,4 +1780,8 @@
 &                                td2n(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","2 m dew point temperature for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
 
             ncret = nf90_def_var(ncid, 'lsprecn', NF90_FLOAT, &
@@ -1412,4 +1796,7 @@
 &                                lsprecn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","large scale total precipitation for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","mm h**-1")
 
             ncret = nf90_def_var(ncid, 'convprecn', NF90_FLOAT, &
@@ -1424,4 +1811,7 @@
 &                                convprecn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","convective precipitation for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","mm h**-1")
 
             ncret = nf90_def_var(ncid, 'sshfn', NF90_FLOAT, &
@@ -1436,4 +1826,7 @@
 &                                sshfn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","surface sensible heat flux for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","J m**-2")
 
             ncret = nf90_def_var(ncid, 'ssrn', NF90_FLOAT, &
@@ -1448,4 +1841,7 @@
 &                                ssrn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","surface solar radiation for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","J m**-2")
 
             ncret = nf90_def_var(ncid, 'surfstrn', NF90_FLOAT, &
@@ -1460,4 +1856,8 @@
 &                                surfstrn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","surface stress for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","N m**-2 s")
+
 
             ncret = nf90_def_var(ncid, 'ustarn', NF90_FLOAT, &
@@ -1472,4 +1872,7 @@
 &                                ustarn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","friction velocity for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
 
             ncret = nf90_def_var(ncid, 'wstarn', NF90_FLOAT, &
@@ -1484,4 +1887,8 @@
 &                                wstarn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","convective velocity scale for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-1")
+
 
             ncret = nf90_def_var(ncid, 'hmixn', NF90_FLOAT, &
@@ -1496,4 +1903,7 @@
 &                                hmixn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","mixing height for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m")
 
             ncret = nf90_def_var(ncid, 'tropopausen', NF90_FLOAT, &
@@ -1508,4 +1918,7 @@
 &                                tropopausen(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","altitude of thermal tropopause for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m")
 
             ncret = nf90_def_var(ncid, 'olin', NF90_FLOAT, &
@@ -1520,4 +1933,7 @@
 &                                olin(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","inverse Obukhov length (1/L) for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m")
 
             ncret = nf90_def_var(ncid, 'diffkn', NF90_FLOAT, &
@@ -1532,4 +1948,10 @@
 &                                diffkn(0:nxmaxn-1, 0:nymaxn-1, 1, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","diffusion coefficient at reference height &
+                                 for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m**2 s**-1")
+
+
 
             dim4dids = (/nxmaxn_dimid, nymaxn_dimid, maxspec_dimid, maxnests_dimid/)
@@ -1548,5 +1970,7 @@
 &                                vdepn(0:nxmaxn-1, 0:nymaxn-1, 1:maxspec, cm_index, 1:maxnests))
             call handle_nf90_err(ncret)
-
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","deposition velocity for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","m s**-2")
 
 !            PRINT *, 'SUM(psn): ', SUM(psn(:,:,:,cm_index,:))
@@ -1570,4 +1994,9 @@
 &                                xresoln(0:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","factor by which the resolutions in the nests &
+                                  is enhanced compared to mother grid")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'yresoln', NF90_FLOAT, &
@@ -1582,4 +2011,9 @@
 &                                yresoln(0:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","factor by which the resolutions in the nests &
+                                  is enhanced compared to mother grid")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             dim1dids = (/maxnests_dimid/)
@@ -1596,4 +2030,9 @@
 &                                xln(1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","lower left corner x  point of nested grids &
+                                in grid coordinates of mother grid")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'yln', NF90_FLOAT, &
@@ -1608,4 +2047,9 @@
 &                                yln(1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","lower left corner y  point of nested grids &
+                                in grid coordinates of mother grid")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'xrn', NF90_FLOAT, &
@@ -1620,4 +2064,9 @@
 &                                xrn(1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","upper right corner x  point of nested grids &
+                                in grid coordinates of mother grid")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'yrn', NF90_FLOAT, &
@@ -1632,4 +2081,8 @@
 &                                yrn(1:maxnests))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","upper right corner y  point of nested grids &
+                                in grid coordinates of mother grid")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
 !            PRINT *, 'SUM(yresoln): ', SUM(yresoln)
@@ -1649,4 +2102,8 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                southpolemap(:))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable to define stereographic projections")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'northpolemap', NF90_FLOAT, &
@@ -1658,5 +2115,7 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                northpolemap(:))
-
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable to define stereographic projections")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
 
@@ -1667,4 +2126,8 @@
             ncret = nf90_put_var(ncid, ncvarid, logical2integer(xglobal))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","logical variable  T for global fields, F for limited area fields")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'sglobal', NF90_INT, ncvarid)
@@ -1672,4 +2135,8 @@
             ncret = nf90_put_var(ncid, ncvarid, logical2integer(sglobal))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","logical variable  T if domain extends towards south pole")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'nglobal', NF90_INT, ncvarid)
@@ -1677,4 +2144,8 @@
             ncret = nf90_put_var(ncid, ncvarid, logical2integer(nglobal))
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","logical variable T if domain extends towards north pole")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'switchnorthg', NF90_FLOAT, ncvarid)
@@ -1682,4 +2153,8 @@
             ncret = nf90_put_var(ncid, ncvarid, switchnorthg)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","degrees for use polar stereographic grid north of switchnorth")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'switchsouthg', NF90_FLOAT, ncvarid)
@@ -1687,4 +2162,7 @@
             ncret = nf90_put_var(ncid, ncvarid, switchsouthg)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","degrees for use polar stereographic grid south of switchshouth")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
 
@@ -1707,4 +2185,8 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                pconv(:))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","pressure at the input model levels for the convection calculation")
+            ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
 
             ncret = nf90_def_var(ncid, 'dpr', NF90_FLOAT, &
@@ -1716,4 +2198,8 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                dpr(:))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","pressure difference around tconv(k)")
+            ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
 
             ncret = nf90_def_var(ncid, 'pconv_hpa', NF90_FLOAT, &
@@ -1725,4 +2211,9 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                pconv_hpa(:))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","pressure at the input model levels &
+                                 in hPa for the convection calculation")
+            ncret = nf90_put_att(ncid, ncvarid, "units","hPa")
+
 
             ncret = nf90_def_var(ncid, 'ft', NF90_FLOAT, &
@@ -1734,4 +2225,8 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                ft(:))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable for the convection scheme")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'fq', NF90_FLOAT, &
@@ -1743,4 +2238,8 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                fq(:))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable for the convection scheme")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
 
             ncret = nf90_def_var(ncid, 'sub', NF90_FLOAT, &
@@ -1753,4 +2252,9 @@
 &                                sub(:))
 
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","auxiliary variable for the convection scheme")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
+
+
             dim1dids = (/na_dimid/)
 
@@ -1763,4 +2267,8 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                phconv(:))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","pressure between levels k-1 and k")
+            ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
 
             ncret = nf90_def_var(ncid, 'phconv_hpa', NF90_FLOAT, &
@@ -1772,4 +2280,7 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                phconv_hpa(:))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","pressure between levels k-1 and k in hPa")
+            ncret = nf90_put_att(ncid, ncvarid, "units","hPa")
 
             ncret = nf90_def_var(ncid, 'tconv', NF90_FLOAT, &
@@ -1781,4 +2292,9 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                tconv(:))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","temperature at the input model levels &
+                                for the convection calculation")
+            ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
 
             ncret = nf90_def_var(ncid, 'qconv', NF90_FLOAT, &
@@ -1790,4 +2306,9 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                qconv(:))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","specific humidity at the input model  &
+                                levels for the convection calculation")
+            ncret = nf90_put_att(ncid, ncvarid, "units","none")
+
 
             ncret = nf90_def_var(ncid, 'qsconv', NF90_FLOAT, &
@@ -1799,4 +2320,9 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                qsconv(:))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","Saturation water vapor specific humidity &
+                                for the convection calculation")
+            ncret = nf90_put_att(ncid, ncvarid, "units","kg/kg")
+
 
             ! New dimensions
@@ -1811,4 +2337,8 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                fmass(:,:))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","mass in the grid column for the convection calculation")
+            ncret = nf90_put_att(ncid, ncvarid, "units","kg")
+
 
             ncret = nf90_def_var(ncid, 'fmassfrac', NF90_FLOAT, &
@@ -1820,4 +2350,9 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                fmassfrac(:,:))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","mass fraction due to  the convection &
+                                displacement from level k to level k")
+            ncret = nf90_put_att(ncid, ncvarid, "units","kg")
+
 
 
@@ -1833,4 +2368,7 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                cbaseflux(0:nxmax-1,0:nymax-1))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","cloudbase massflux due to the convection")
+            ncret = nf90_put_att(ncid, ncvarid, "units","(kg/m**2)/s")
 
             ! New dimensions
@@ -1845,4 +2383,7 @@
             ncret = nf90_put_var(ncid, ncvarid, &
 &                                cbasefluxn(0:nxmaxn-1,0:nymaxn-1,1:maxnests))
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","cloudbase massflux due to the convection for the nested domain")
+            ncret = nf90_put_att(ncid, ncvarid, "units","(kg/m**2)/s")
 
 
@@ -1852,4 +2393,8 @@
             ncret = nf90_put_var(ncid, ncvarid, psconv)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","surface pressure for the convection calculation")
+            ncret = nf90_put_att(ncid, ncvarid, "units","Pa")
+
 
             ncret = nf90_def_var(ncid, 'tt2conv', NF90_FLOAT, ncvarid)
@@ -1857,4 +2402,8 @@
             ncret = nf90_put_var(ncid, ncvarid, tt2conv)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","2 m temperature for the convection calculation")
+            ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
 
             ncret = nf90_def_var(ncid, 'td2conv', NF90_FLOAT, ncvarid)
@@ -1862,4 +2411,8 @@
             ncret = nf90_put_var(ncid, ncvarid, td2conv)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","2 m dew point temperature for the convection calculation")
+            ncret = nf90_put_att(ncid, ncvarid, "units","K")
+
 
             ncret = nf90_def_var(ncid, 'nconvlev', NF90_INT, ncvarid)
@@ -1867,4 +2420,7 @@
             ncret = nf90_put_var(ncid, ncvarid, nconvlev)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","number of levels for convection")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
             ncret = nf90_def_var(ncid, 'nconvtop', NF90_INT, ncvarid)
@@ -1872,4 +2428,7 @@
             ncret = nf90_put_var(ncid, ncvarid, nconvtop)
             call handle_nf90_err(ncret)
+! attributes
+            ncret = nf90_put_att(ncid, ncvarid, "description","upper levels for convection")
+            ncret = nf90_put_att(ncid, ncvarid, "units"," ")
 
 !            PRINT *, 'SUM(pconv): ', SUM(pconv)
Index: flexpart_code/grib2nc4/Makefile
===================================================================
--- flexpart_code/grib2nc4/Makefile	(revision 87d96843f7f85bb599c0f7cdfb7f1f42731b9eaf)
+++ flexpart_code/grib2nc4/Makefile	(revision 8624a75e20c0cc6dff42b4eb97182a138f9aa1e7)
@@ -11,19 +11,23 @@
 NETCDFF = /opt/netcdf-fortran-4.4.3
 NETCDF = /opt/netcdf-c-4.4.0
-
+#GRIBAPI = /usr/local/grib-api
+#HDF5 = /usr/local/hdf5-1.8.16
+#NETCDFF = /usr/local/netcdf-fortran-4.4.3
+#NETCDF = /usr/local/netcdf-c-4.4.0
 
 
 
 BINARY = grib2nc4
-OBJS = processmetfields.o
-FPMODOBJS = par_mod.o com_mod.o class_vtable_mod.o cmapf_mod.o conv_mod.o 
-FLXPRTOBJS = detectformat.o grib2check.o shift_field_0.o gridcheck.o \
-	     readwind.o readwind_nests.o calcpar.o calcpar_nests.o \
-             shift_field.o pbl_profile.o scalev.o obukhov.o \
-             richardson.o ew.o getvdep.o calcpv.o obukhov_gfs.o \
-             richardson_gfs.o getvdep_nests.o calcpv_nests.o psim.o psih.o \
-             qvsat.o caldate.o getrb.o raerod.o getrc.o partdep.o \
-             verttransform.o verttransform_nests.o readwind_gfs.o \
-             calcpar_gfs.o verttransform_gfs.o gridcheck_gfs.o
+CMP_BINARY = nc4cmp
+OBJS = processmetfields.o verttransform_grib2nc4_ecmwf.o verttransform_grib2nc4_gfs.o
+FPMODOBJS = ${FLEXPART_SRC}/par_mod.o ${FLEXPART_SRC}/com_mod.o ${FLEXPART_SRC}/class_vtable_mod.o ${FLEXPART_SRC}/cmapf_mod.o ${FLEXPART_SRC}/conv_mod.o 
+FLXPRTOBJS = ${FLEXPART_SRC}/detectformat.o ${FLEXPART_SRC}/grib2check.o ${FLEXPART_SRC}/shift_field_0.o ${FLEXPART_SRC}/gridcheck.o \
+	     ${FLEXPART_SRC}/readwind.o ${FLEXPART_SRC}/readwind_nests.o ${FLEXPART_SRC}/calcpar.o ${FLEXPART_SRC}/calcpar_nests.o \
+             ${FLEXPART_SRC}/shift_field.o ${FLEXPART_SRC}/pbl_profile.o ${FLEXPART_SRC}/scalev.o ${FLEXPART_SRC}/obukhov.o \
+             ${FLEXPART_SRC}/richardson.o ${FLEXPART_SRC}/ew.o ${FLEXPART_SRC}/getvdep.o ${FLEXPART_SRC}/calcpv.o ${FLEXPART_SRC}/obukhov_gfs.o \
+             ${FLEXPART_SRC}/richardson_gfs.o ${FLEXPART_SRC}/getvdep_nests.o ${FLEXPART_SRC}/calcpv_nests.o ${FLEXPART_SRC}/psim.o ${FLEXPART_SRC}/psih.o \
+             ${FLEXPART_SRC}/qvsat.o ${FLEXPART_SRC}/caldate.o ${FLEXPART_SRC}/getrb.o ${FLEXPART_SRC}/raerod.o ${FLEXPART_SRC}/getrc.o ${FLEXPART_SRC}/partdep.o \
+             ${FLEXPART_SRC}/verttransform.o ${FLEXPART_SRC}/verttransform_nests.o ${FLEXPART_SRC}/readwind_gfs.o \
+             ${FLEXPART_SRC}/calcpar_gfs.o ${FLEXPART_SRC}/verttransform_gfs.o ${FLEXPART_SRC}/gridcheck_gfs.o
 
 
@@ -32,14 +36,14 @@
 
 INCLUDES_NETCDF = -I${NETCDFF}/include
-INCLUDES = -I${GRIBAPI}/include ${INCLUDES_NETCDF}
+INCLUDES = -I${GRIBAPI}/include ${INCLUDES_NETCDF}  -I${FLEXPART_SRC}
 
 
 ### NetCDF link flags - use the first one for dynamic libs, the second
 ### one for static libs
-#LDFLAGS_NETCDF = -L${NETCDFF}/lib -lnetcdff -L${NETCDF}/lib -lnetcdf 
-LDFLAGS_NETCDF=-static -L${NETCDFF}/lib -lnetcdff -L${NETCDF}/lib -lnetcdf -lnetcdf -L${HDF5}/lib -lhdf5_fortran -lhdf5_hl -lhdf5hl_fortran -lhdf5 -ldl -lz
+LDFLAGS_NETCDF = -L${NETCDFF}/lib -lnetcdff -L${NETCDF}/lib -lnetcdf 
+#LDFLAGS_NETCDF=-static -L${NETCDFF}/lib -lnetcdff -L${NETCDF}/lib -lnetcdf -lnetcdf -L${HDF5}/lib -lhdf5_fortran -lhdf5_hl -lhdf5hl_fortran -lhdf5 -ldl -lz
 
 
-LDFLAGS = -L${GRIBAPI}/lib -lgrib_api_f90 -lgrib_api ${LDFLAGS_NETCDF} -ljasper
+LDFLAGS = -L${GRIBAPI}/lib -lgrib_api_f90 -lgrib_api ${LDFLAGS_NETCDF} -ljasper -L${FLEXPART_SRC}
 
 
@@ -51,4 +55,11 @@
 ${BINARY}.o : ${BINARY}.F90 fp2nc4io_mod.mod ${FPMODOBJS} Makefile
 	${FC} -c ${BINARY}.F90 ${FFLAGS} ${INCLUDES}
+
+#-----------  NC4 compare ------------------------
+${CMP_BINARY} : ${CMP_BINARY}.o
+	${FC} -o ${CMP_BINARY} ${CMP_BINARY}.o ${LDFLAGS}
+
+${CMP_BINARY}.o : ${CMP_BINARY}.F90 Makefile
+	${FC} -c ${CMP_BINARY}.F90 ${FFLAGS} ${INCLUDES}
 
 fp2nc4io_mod.mod : ${FPMODOBJS} 
Index: flexpart_code/grib2nc4/README
===================================================================
--- flexpart_code/grib2nc4/README	(revision 496c6079149f12fbb3e7ecf1f351470794679655)
+++ flexpart_code/grib2nc4/README	(revision 8624a75e20c0cc6dff42b4eb97182a138f9aa1e7)
@@ -7,5 +7,5 @@
          M. Harustak
 
-Last Update: 29 May 2016
+Last Update: 05 July 2017
 
 ======================================================================
@@ -19,4 +19,8 @@
 data for computations.  It then takes the data which has been stored in
 FLEXPART global arrays and selectively writes it to a NetCDF4 file.
+
+Optional parameter allows user to request arrays calculated for specified
+location instead of default algorithm selecting lowest left corner of meteorological
+domain with a surface pressure above 100000Pa.
 
 This distribution is integrated into the FLEXPART source code, in a grib2nc4
@@ -42,4 +46,12 @@
 is good, but only whether data was stored correctly in the NC4 files.
 
+To be able to quickly compare variables in NetCDF files, simple utility
+nc4cmp was introduced. It compares specified variable in two NetCDF files
+(with dimension up to 3) and it prints average and maximal difference between
+the values. It allows for specification of optional tolerance (in %).
+Difference is calculated as abs((value1-value2)/value1)*100.
+It is calculated for all values of variable multidimensional array.
+
+
 ======================================================================
 2. Quick-Start
@@ -64,4 +76,6 @@
 Simply type "make" (noting that you may see some warnings related to
 variables for nests) to produce grib2nc4.
+
+NetCDF compare utility nc4cmp can be built using command "make nc4cmp".
 
 Before trying to run, you might want to try "make test" to make sure that
@@ -81,4 +95,9 @@
 as the second argument.
 
+Usage:
+
+  ./grib2nc4 <inpath> <outpath> [lat=value lon=value] [optional varnames]
+
+
 Sample usage:
 
@@ -92,4 +111,27 @@
 
   ./grib2nc4 GD15051200 def9_nc1p0.nc4  q  w
+
+Calculation at user specified location can be requested by adding coord parameter
+to the command line:
+
+  ./grib2nc4 GD15051200 def9_nc1p0.nc4 lon=20 lat=45
+
+To compare two NetCDF files, nc4cmp can be used.
+
+Usage:
+
+  ./nc4cmp <file1> <file2> <variable> [optional tolerance in %]
+
+Sample usage:
+
+  ./nc4cmp def9_nc1p0.nc4 def9_nc1p0_reference.nc4 U
+
+Sample usage with tolerance 1%:
+
+  ./nc4cmp def9_nc1p0.nc4 def9_nc1p0_reference.nc4 U 1.0
+
+This will extract variable U from both NetCDF files, print and compare dimensions,
+compare values and print average and maximal difference.
+
 
 ======================================================================
Index: flexpart_code/grib2nc4/fp2nc4io_mod.F90
===================================================================
--- flexpart_code/grib2nc4/fp2nc4io_mod.F90	(revision 87d96843f7f85bb599c0f7cdfb7f1f42731b9eaf)
+++ flexpart_code/grib2nc4/fp2nc4io_mod.F90	(revision 8624a75e20c0cc6dff42b4eb97182a138f9aa1e7)
@@ -157,4 +157,7 @@
         ncfunc_retval = nf90_def_var(ncid, 'height', NF90_DOUBLE, &
 &                                    z_dimid, varid)
+! attributes
+            ncfunc_retval = nf90_put_att(ncid, varid, "description","height of the FLEXPART model levels")
+            ncfunc_retval = nf90_put_att(ncid, varid, "units","m a.g.l")
 
         ncfunc_retval = nf90_def_var_deflate(ncid, varid,   &
@@ -169,13 +172,29 @@
         ncfunc_retval = nf90_def_var(ncid, 'dx', NF90_DOUBLE, varid)
         ncfunc_retval = nf90_put_var(ncid, varid, dx)
+! attributes
+            ncfunc_retval = nf90_put_att(ncid, varid, "description","grid distance in x direction")
+            ncfunc_retval = nf90_put_att(ncid, varid, "units","degrees")
+
 
         ncfunc_retval = nf90_def_var(ncid, 'dy', NF90_DOUBLE, varid)
         ncfunc_retval = nf90_put_var(ncid, varid, dy)
+! attributes
+            ncfunc_retval = nf90_put_att(ncid, varid, "description","grid distance in y direction")
+            ncfunc_retval = nf90_put_att(ncid, varid, "units","degrees")
+
 
         ncfunc_retval = nf90_def_var(ncid, 'xlon0', NF90_DOUBLE, varid)
         ncfunc_retval = nf90_put_var(ncid, varid, xlon0)
+! attributes
+            ncfunc_retval = nf90_put_att(ncid, varid, "description","longitude of the lowest left corner")
+            ncfunc_retval = nf90_put_att(ncid, varid, "units","degrees")
+
 
         ncfunc_retval = nf90_def_var(ncid, 'ylat0', NF90_DOUBLE, varid)
         ncfunc_retval = nf90_put_var(ncid, varid, ylat0)
+! attributes
+            ncfunc_retval = nf90_put_att(ncid, varid, "description","latitude of the lowest left corner")
+            ncfunc_retval = nf90_put_att(ncid, varid, "units","degrees")
+
 
         ! All done, close the NetCDF file
@@ -296,16 +315,40 @@
             ncfunc_retval = nf90_put_var(ncid, varid, &
 &                                        uu(0:nx-1, 0:ny-1, 1:nz, 1))
+! attributes
+            ncfunc_retval = nf90_put_att(ncid, varid, "description","U component of wind in the X[horizontal] direction")
+            ncfunc_retval = nf90_put_att(ncid, varid, "units","m s**-1")
+
+
         ELSEIF (nc_varname == 'V') THEN
             ncfunc_retval = nf90_put_var(ncid, varid, &
 &                                        vv(0:nx-1, 0:ny-1, 1:nz, 1))
+! attributes
+            ncfunc_retval = nf90_put_att(ncid, varid, "description","V component of wind in the Y[horizontal] direction")
+            ncfunc_retval = nf90_put_att(ncid, varid, "units","m s**-1")
+
+
         ELSEIF (nc_varname == 'T') THEN
             ncfunc_retval = nf90_put_var(ncid, varid, &
 &                                        tt(0:nx-1, 0:ny-1, 1:nz, 1))
+! attributes
+            ncfunc_retval = nf90_put_att(ncid, varid, "description","temperature")
+            ncfunc_retval = nf90_put_att(ncid, varid, "units","k")
+
         ELSEIF (nc_varname == 'W') THEN
             ncfunc_retval = nf90_put_var(ncid, varid, &
 &                                        ww(0:nx-1, 0:ny-1, 1:nz, 1))
+! attributes
+            ncfunc_retval = nf90_put_att(ncid, varid, "description","wind component in the Z[vertical] direction")
+            ncfunc_retval = nf90_put_att(ncid, varid, "units","m s**-1")
+
+
         ELSEIF (nc_varname == 'Q') THEN
             ncfunc_retval = nf90_put_var(ncid, varid, &
 &                                        qv(0:nx-1, 0:ny-1, 1:nz, 1))
+! attributes
+            ncfunc_retval = nf90_put_att(ncid, varid, "description","specific humidity")
+            ncfunc_retval = nf90_put_att(ncid, varid, "units"," ")
+
+
         ELSE
             PRINT *,
Index: flexpart_code/grib2nc4/grib2nc4.F90
===================================================================
--- flexpart_code/grib2nc4/grib2nc4.F90	(revision 496c6079149f12fbb3e7ecf1f351470794679655)
+++ flexpart_code/grib2nc4/grib2nc4.F90	(revision 8624a75e20c0cc6dff42b4eb97182a138f9aa1e7)
@@ -11,4 +11,9 @@
     !        May 2016                                                        *
     !*************************************************************************
+   !   M. Harustak                                                          *
+   !   -) modification to generate the output in single precission          *
+   !   -) possibility to add a lat lon selection to obtain the met variables*
+   !      in the vertical levels defined in that location                   *
+   !*************************************************************************
 
     USE par_mod
@@ -21,11 +26,12 @@
     IMPLICIT NONE
 
-    LOGICAL :: metfile_exists   
+    LOGICAL :: metfile_exists, coordinates_provided, lat_provided, lon_provided
     INTEGER :: i, j, k
     INTEGER :: num_optional_vars, num_vars
     INTEGER, PARAMETER :: DEFLATE_LEVEL = 2  ! Compression level (0-9)
-    CHARACTER(LEN=512) :: met_filepath, netcdf4_filepath
+    CHARACTER(LEN=512) :: met_filepath, netcdf4_filepath, param_str, coord_name_str, coord_val_str
     CHARACTER, DIMENSION(:), ALLOCATABLE :: vars_list  ! list of variables
-    
+    INTEGER :: coordX, coordY, stat
+    REAL :: coord_lat, coord_lon
     INTEGER :: metdata_format = UNKNOWN_METDATA  ! From FP par_mod
 
@@ -34,5 +40,5 @@
     ! Read in mandatory arguments
     IF (IARGC() < 2) THEN
-        PRINT *, 'Usage: grib2netcdf4 <inpath> <outpath> [optional varnames]'
+        PRINT *, 'Usage: grib2netcdf4 <inpath> <outpath> [lon=X lat=Y] [optional varnames]'
         STOP
     ELSE
@@ -48,21 +54,45 @@
     ! First, get the number of optional args and allocate vars_list,
     ! and fill the first three elements
-    IF (IARGC() > 2) THEN
-        num_optional_vars = IARGC() - 2
-    ELSE
-        num_optional_vars = 0
-    ENDIF
+    coordinates_provided = .FALSE.
+    lat_provided = .FALSE.
+    lon_provided = .FALSE.
+    ALLOCATE( vars_list(IARGC()+3),stat=stat )
 
-    num_vars = num_optional_vars + 3
-    ALLOCATE( vars_list(num_vars) )
     vars_list(1) = 'u'
     vars_list(2) = 'v'
     vars_list(3) = 't'
         
-    ! Read in optional variable arguments, starting at element 4 of vars_list
-    IF (IARGC() > 2) THEN
-        DO i=1,num_optional_vars
-            CALL GETARG( i+2, vars_list(i+3) )
-        ENDDO 
+    num_vars = 3
+    DO i=3,IARGC()
+        CALL GETARG(i,param_str)
+        param_str = TRIM(param_str)
+        j = SCAN(param_str,"=")
+        if (j>1) then
+            coord_name_str=param_str(1:j-1)
+            coord_val_str=param_str(j+1:)
+            IF ( coord_name_str == "lat" .or. coord_name_str == "LAT" ) THEN
+                read(coord_val_str,*,iostat=stat) coord_lat
+                if ( stat == 0 ) then
+                    lat_provided = .TRUE.
+                else
+                    print *, "Incorrect coordinates: ", coord_val_str
+                    stop
+                endif
+            ELSE IF ( coord_name_str == "lon" .or. coord_name_str == "LON" ) THEN
+                read(coord_val_str,*,iostat=stat) coord_lon
+                if ( stat == 0 ) then
+                    lon_provided = .TRUE.
+                else
+                    print *, "Incorrect coordinates: ", coord_val_str
+                    stop
+                endif
+            ENDIF
+        else
+            num_vars = num_vars + 1
+            vars_list(num_vars) = param_str
+        endif
+    ENDDO
+    IF (lat_provided .AND. lon_provided) THEN
+        coordinates_provided = .TRUE.
     ENDIF
 
@@ -126,6 +156,16 @@
     !CALL gridcheck_nests
 
+    ! If the coordinates are provided, then we need to obtain
+    !  the corresponding grid indexes since this is what verttransform needs
+    if ( coordinates_provided ) then
+        coordX = (coord_lon-xlon0)/dx
+        coordY = (coord_lat-ylat0)/dy
+        print *, "Coordinates: "
+        print *, "lon: ", coord_lon, ", lat: ",coord_lat
+        print *, "x: ", coordX, ", y: ", coordY
+    endif
+
     PRINT *, 'Calling processmetfields()...'
-    call processmetfields( 1, metdata_format)
+    call processmetfields( 1, metdata_format, coordinates_provided, coordX, coordY)
 
     PRINT *, 'Calling fp2nc4io_dump()...'
Index: flexpart_code/grib2nc4/nc4cmp.F90
===================================================================
--- flexpart_code/grib2nc4/nc4cmp.F90	(revision 8624a75e20c0cc6dff42b4eb97182a138f9aa1e7)
+++ flexpart_code/grib2nc4/nc4cmp.F90	(revision 8624a75e20c0cc6dff42b4eb97182a138f9aa1e7)
@@ -0,0 +1,244 @@
+PROGRAM nc4cmp
+
+    !*************************************************************************
+    !  This program compares variable in two NetCDF files                    *
+    ! (with dimension up to 3) and it prints average and maximal difference  *
+    ! between the values. It allows for specification of optional tolerance  *
+    !                                                                        *
+    ! Usage:                                                                 *
+    !                                                                        *
+    !     ./nc4cmp <file1> <file2> <variable> [optional tolerance in %]      *
+    !                                                                        *
+    !        M. Harustak                                                     *
+    !                                                                        *
+    !        March 2017                                                      *
+    !*************************************************************************
+
+
+    USE netcdf
+
+    IMPLICIT NONE
+
+    LOGICAL :: file1_exists, file2_exists
+    CHARACTER(LEN=512) :: file1, file2, precision_str, variable
+    CHARACTER(LEN=NF90_MAX_NAME) :: returned_name
+    REAL :: tolerance
+    INTEGER :: ncid, varid, retval, xtype, numDims1, numDims2, i, j, k, stat
+    INTEGER, dimension(NF90_MAX_VAR_DIMS) :: varDimIDs1, varDimIDs2, dims1, dims2
+    double precision :: variable_array_0d_1, variable_array_0d_2
+    double precision, allocatable, dimension(:) :: variable_array_1d_1, variable_array_1d_2 
+    double precision, allocatable, dimension(:,:) :: variable_array_2d_1, variable_array_2d_2
+    double precision, allocatable, dimension(:,:,:) :: variable_array_3d_1, variable_array_3d_2
+    double precision :: v1, v2
+    double precision :: sum_diff, max_diff
+    integer :: count_diff
+    logical :: different = .FALSE.
+
+    ! Read in mandatory arguments
+    IF (IARGC() < 3) THEN
+        PRINT *, 'Usage: nc4cmp <file1> <file2> <variable> [optional tolerance in %]'
+        STOP
+    ELSE
+        CALL GETARG(1, file1) 
+        CALL GETARG(2, file2) 
+        CALL GETARG(3, variable)
+    ENDIF
+
+    IF (IARGC() > 3) THEN
+        CALL GETARG( 4, precision_str)
+        READ(precision_str, '(f10.5)') tolerance
+    ELSE
+        tolerance = 0
+    ENDIF
+    print *,"Tolerance ",tolerance,"%"
+
+    INQUIRE( FILE=file1, EXIST=file1_exists )
+    IF ( file1_exists ) THEN
+        PRINT *, 'File 1: ', TRIM(file1)
+    ELSE
+        PRINT *, 'Unable to find file: ', TRIM(file1)
+        STOP
+    ENDIF
+    INQUIRE( FILE=file2, EXIST=file2_exists )
+    IF ( file2_exists ) THEN
+        PRINT *, 'File 2: ', TRIM(file2)
+    ELSE
+        PRINT *, 'Unable to find file: ', TRIM(file2)
+        STOP
+    ENDIF
+
+    print *, "Opening ",TRIM(file1)," for reading"
+    retval = nf90_open(file1, NF90_NOWRITE, ncid)
+    retval = nf90_inq_varid(ncid, TRIM(variable) , varid)
+    retval = nf90_inquire_variable(ncid, varid, ndims = numDims1, dimids = varDimIDs1)
+    print *, "Reading variable ",TRIM(variable), ", no of dimensions: ",numDims1
+    do i=1,numDims1
+        retval = nf90_inquire_dimension(ncid, varDimIDs1(i),len=j)
+        dims1(i) = j
+        print *,"  - dimension ",i,": ", j
+    enddo
+    if ( numDims1 == 0) then
+        retval = nf90_get_var(ncid, varid, variable_array_0d_1)
+    else if ( numDims1 == 1) then
+        allocate(variable_array_1d_1(dims1(1)),stat=stat)
+        retval = nf90_get_var(ncid, varid, variable_array_1d_1)
+    else if ( numDims1 == 2) then
+        allocate(variable_array_2d_1(dims1(1),dims1(2)),stat=stat)
+        retval = nf90_get_var(ncid, varid, variable_array_2d_1)
+    else if ( numDims1 == 3) then
+        allocate(variable_array_3d_1(dims1(1),dims1(2),dims1(3)),stat=stat)
+        retval = nf90_get_var(ncid, varid, variable_array_3d_1)
+    endif
+    retval = nf90_close(ncid)
+    print *, "...done"
+
+    print *, "Opening ",TRIM(file2)," for reading"
+    retval = nf90_open(file2, NF90_NOWRITE, ncid)
+    retval = nf90_inq_varid(ncid, TRIM(variable) , varid)
+    retval = nf90_inquire_variable(ncid, varid, ndims = numDims2, dimids = varDimIDs2)
+    print *, "Reading variable ",TRIM(variable), ", no of dimensions: ",numDims2
+    do i=1,numDims2
+        retval = nf90_inquire_dimension(ncid, varDimIDs2(i),len=j)
+        dims2(i) = j 
+        print *,"  - dimension ",i,": ", j
+    enddo
+    if ( numDims2 == 0) then
+        retval = nf90_get_var(ncid, varid, variable_array_0d_2)
+    else if ( numDims2 == 1) then
+        allocate(variable_array_1d_2(dims2(1)),stat=stat)
+        retval = nf90_get_var(ncid, varid, variable_array_1d_2)
+    else if ( numDims2 == 2) then
+        allocate(variable_array_2d_2(dims2(1),dims2(2)),stat=stat)
+        retval = nf90_get_var(ncid, varid, variable_array_2d_2)
+    else if ( numDims2 == 3) then
+        allocate(variable_array_3d_2(dims2(1),dims2(2),dims2(3)),stat=stat)
+        retval = nf90_get_var(ncid, varid, variable_array_3d_2)
+    endif
+    retval = nf90_close(ncid)
+    print *, "...done"
+
+    if (numDims1 /= numDims2) THEN
+        print *, "Number of dimensions differs..."
+        STOP
+    endif
+    do i=1,numDims1
+        if (dims1(i) /= dims2(i)) then
+            print *, "Dimension ",i," differs..."
+            STOP
+        endif
+    enddo
+
+    if (numDims1>3) then
+        print *,"Only up to 3 dimensions supported..." 
+        STOP
+    endif
+
+    sum_diff = 0.0
+    count_diff = 0
+
+    if (numDims1==0) then
+        v1 = variable_array_0d_1
+        v2 = variable_array_0d_2
+        print *, "Value 1: ",v1
+        print *, "Value 2: ",v2
+        if ( ( v1 == 0 .AND. v2 /= 0 ) .OR. ( v1 /= 0 .AND. v2 == 0 ) ) then
+            print *," === variables differ ==="
+            STOP
+        else if ( v1 == 0 .AND. v2 == 0) then
+            print *," === variables are within tolerance ==="
+            STOP
+        else if ( abs((v1-v2)/v1*100) > tolerance ) then 
+            print *," === variables differ ==="
+            STOP
+        else 
+            print *," === variables are within tolerance ==="
+            STOP
+        endif
+    endif
+    if (numDims1==1) then
+        do i=1,dims1(1)
+            v1 = variable_array_1d_1(i)
+            v2 = variable_array_1d_2(i)
+            if ( ( v1 == 0 .AND. v2 /= 0 ) .OR. ( v1 /= 0 .AND. v2 == 0 ) ) then
+                different = .TRUE. 
+            else if (v1/=0 .AND. v2/=0 .AND. abs((v1-v2)/v1*100) > tolerance ) then
+                different = .TRUE.
+            endif
+            if ( v1 /= 0) then 
+                sum_diff = sum_diff+abs((v1-v2)/v1*100)
+                count_diff = count_diff+1
+                if (abs((v1-v2)/v1*100)>max_diff) then
+                    max_diff=abs((v1-v2)/v1*100)
+                endif
+            endif
+        enddo
+        if ( different ) then
+            print *," === variables differ ==="
+        else
+            print *," === variables are within tolerance ==="
+        endif
+        print *," === avg(diff%): ", sum_diff/count_diff
+        print *," === max(diff%): ",max_diff
+        STOP
+    endif
+    if (numDims1==2) then
+        do i=1,dims1(1)
+            do j=1,dims1(2)
+                v1 = variable_array_2d_1(i,j)
+                v2 = variable_array_2d_2(i,j)
+                if ( ( v1 == 0 .AND. v2 /= 0 ) .OR. ( v1 /= 0 .AND. v2 == 0 ) ) then
+                    different = .TRUE.
+                else if (v1/=0 .AND. v2/=0 .AND. abs((v1-v2)/v1*100) > tolerance ) then
+                    different = .TRUE.
+                endif
+                if (v1/=0) then
+                    sum_diff = sum_diff+abs((v1-v2)/v1*100)
+                    count_diff = count_diff+1
+                    if (abs((v1-v2)/v1*100)>max_diff) then
+                        max_diff=abs((v1-v2)/v1*100)
+                    endif
+                endif
+            enddo
+        enddo
+        if ( different ) then
+            print *," === variables differ ==="
+        else
+            print *," === variables are within tolerance ==="
+        endif
+        print *," === avg(diff%): ", sum_diff/count_diff
+        print *," === max(diff%): ",max_diff
+        STOP
+    endif
+    if (numDims1==3) then
+        do i=1,dims1(1)
+            do j=1,dims1(2)
+                do k=1,dims1(3)
+                    v1 = variable_array_3d_1(i,j,k)
+                    v2 = variable_array_3d_2(i,j,k)
+                    if ( ( v1 == 0 .AND. v2 /= 0 ) .OR. ( v1 /= 0 .AND. v2 == 0 ) ) then
+                        different = .TRUE.
+                    else if (v1/=0 .AND. v2/=0 .AND. abs((v1-v2)/v1*100) > tolerance ) then
+                        different = .TRUE.
+                    endif
+                    if (v1/=0) then
+                        sum_diff = sum_diff+abs((v1-v2)/v1*100)
+                        count_diff = count_diff+1
+                        if (abs((v1-v2)/v1*100)>max_diff) then
+                            max_diff=abs((v1-v2)/v1*100)
+                        endif
+                    endif
+                enddo
+            enddo
+        enddo
+        if ( different ) then
+            print *," === variables differ ==="
+        else
+            print *," === variables are within tolerance ==="
+        endif
+        print *," === avg(diff%): ", sum_diff/count_diff
+        print *," === max(diff%): ",max_diff
+        STOP
+    endif
+
+
+END PROGRAM nc4cmp
Index: flexpart_code/grib2nc4/processmetfields.F90
===================================================================
--- flexpart_code/grib2nc4/processmetfields.F90	(revision 496c6079149f12fbb3e7ecf1f351470794679655)
+++ flexpart_code/grib2nc4/processmetfields.F90	(revision 8624a75e20c0cc6dff42b4eb97182a138f9aa1e7)
@@ -1,3 +1,4 @@
-subroutine processmetfields(ind,metdata_format)
+subroutine processmetfields(ind,metdata_format,coordinates_provided, coordX, coordY)
+
   !                       i     o
   !*****************************************************************************
@@ -37,4 +38,6 @@
   integer :: dumpData
   character(len=512):: fpfname, dumpPath, filename
+  integer :: coordX, coordY
+  logical :: coordinates_provided
 
 
@@ -107,5 +110,9 @@
              call calcpar_ecmwf(memind(1),uuh,vvh,pvh) 
              call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format)
-             call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh)
+             if ( coordinates_provided ) then
+                 call verttransform_grib2nc4_ecmwf(memind(1),uuh,vvh,wwh,pvh,coordX,coordY)
+             else
+                 call verttransform_grib2nc4_ecmwf(memind(1),uuh,vvh,wwh,pvh,-1,-1)
+             endif
              call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
              memtime(1)=wftime(ind)
@@ -116,5 +123,9 @@
              call calcpar_gfs(memind(1),uuh,vvh,pvh)
              call calcpar_nests(memind(1),uuhn,vvhn,pvhn,metdata_format)
-             call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh)
+             if ( coordinates_provided ) then
+                 call verttransform_grib2nc4_gfs(memind(1),uuh,vvh,wwh,pvh,coordX,coordY)
+             else
+                 call verttransform_grib2nc4_gfs(memind(1),uuh,vvh,wwh,pvh,-1,-1)
+             endif
              call verttransform_nests(memind(1),uuhn,vvhn,wwhn,pvhn)
              memtime(1)=wftime(ind)
Index: flexpart_code/grib2nc4/verttransform_grib2nc4_ecmwf.f90
===================================================================
--- flexpart_code/grib2nc4/verttransform_grib2nc4_ecmwf.f90	(revision 8624a75e20c0cc6dff42b4eb97182a138f9aa1e7)
+++ flexpart_code/grib2nc4/verttransform_grib2nc4_ecmwf.f90	(revision 8624a75e20c0cc6dff42b4eb97182a138f9aa1e7)
@@ -0,0 +1,701 @@
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
+!                                                                     *
+! This file is part of FLEXPART.                                      *
+!                                                                     *
+! FLEXPART is free software: you can redistribute it and/or modify    *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or   *
+! (at your option) any later version.                                 *
+!                                                                     *
+! FLEXPART is distributed in the hope that it will be useful,         *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
+! GNU General Public License for more details.                        *
+!                                                                     *
+! You should have received a copy of the GNU General Public License   *
+! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
+!**********************************************************************
+
+subroutine verttransform_grib2nc4_ecmwf(n,uuh,vvh,wwh,pvh,coordx_opt,coordy_opt)
+  !                         i  i   i   i   i
+  !*****************************************************************************
+  !                                                                            *
+  !     This subroutine transforms temperature, dew point temperature and      *
+  !     wind components from eta to meter coordinates.                         *
+  !     The vertical wind component is transformed from Pa/s to m/s using      *
+  !     the conversion factor pinmconv.                                        *
+  !     In addition, this routine calculates vertical density gradients        *
+  !     needed for the parameterization of the turbulent velocities.           *
+  !                                                                            *
+  !     Author: A. Stohl, G. Wotawa                                            *
+  !                                                                            *
+  !     12 August 1996                                                         *
+  !     Update: 16 January 1998                                                *
+  !                                                                            *
+  !     Major update: 17 February 1999                                         *
+  !     by G. Wotawa                                                           *
+  !                                                                            *
+  !     - Vertical levels for u, v and w are put together                      *
+  !     - Slope correction for vertical velocity: Modification of calculation  *
+  !       procedure                                                            *
+  !                                                                            *
+  !*****************************************************************************
+  !  Changes, Bernd C. Krueger, Feb. 2001:
+  !   Variables tth and qvh (on eta coordinates) from common block
+  !*****************************************************************************
+  ! Sabine Eckhardt, March 2007
+  ! added the variable cloud for use with scavenging - descr. in com_mod
+  !*****************************************************************************
+  !  Changes Arnold, D. and Morton, D. (2015):                                 *
+  !   -- description of local and common variables                             *
+  !*****************************************************************************
+  !  July 2017 - M. Harustak                                                   *
+  !  In order to be used by GRIB2NC4, this routine has been modified so that   *
+  !  it has the possibility to get two additional arguments, coordx_opt and    *
+  !  coordy_opt. These two are used to define the location where the FLEXPART  *
+  !  model levels are calculated and used in GRIB2NC4 to obtaine the met data  *
+  !                                                                            *
+  !  The routine works normally if -1 (no value) is passed                     *
+  !*****************************************************************************
+
+
+  use par_mod
+  use com_mod
+  use cmapf_mod, only: cc2gll
+
+  implicit none
+
+  !***********************************************************************
+  ! Subroutine Parameters:                                               *
+  !    input                                                             *
+  ! n                   temporal index for meteorological fields (1 to 3)*  
+  ! uuh,vvh, wwh        wind components in ECMWF model levels            *  
+  ! pvh                 potential vorticity in ECMWF model levels        *
+  integer :: n
+  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
+  integer :: coordx_opt, coordy_opt
+  !***********************************************************************
+
+  !***********************************************************************
+  ! Local variables                                                      *
+  !                                                                      *
+  !  ew                   subroutine/function to calculate saturation    *  
+  !                         water vaport for a given air temperature     *
+  ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition     *
+  ! rain_cloud_above [0,1]  whether there is a raining cloud or not      *  
+  ! kz_inv  inverted indez for kz                                        *
+  ! f_qvsat  Saturation water vapor specific humidity (kg/kg)            *
+  ! pressure                                                             *
+  ! rh       relative humidity                                           *    
+  ! lsp      large scale precip at one grid cell                         *
+  ! convp    convective precip at one grid cell                          *
+  ! uvzlev    height of the eta half-levels                              *
+  ! rhoh     density in model levels                                     *
+  ! pinmconv conversion factor                                           *
+  ! pint     pressure on model levels (using akz, bkz, ps)               *
+  ! tv       virtual temperature                                         *
+  ! tvold,pold      dummy variables to keep previous value               *
+  ! dz1, dz2, dz  differences heights model levels, pressure levels      *    
+  ! ui, vi                                                               *   
+  ! xlon,xlat                                                            *
+  ! xlonr         xlon*pi/180.                                           * 
+  ! dzdx,dzdy     slope corrections                                      *
+  ! dzdx1,dzdx2,dzdy1   slope corrections in each direction              *
+  ! uuaux,vvaux,uupolaux,vvpolaux  auxiliary variables for polar sterero.*
+  ! ddpol,ffpol  for special treatment of poles                          *
+  ! wdummy                                                               *
+  ! wzlev                                                                *
+  ! uvwzlev                                                              *
+
+  integer :: rain_cloud_above,kz_inv
+  real :: f_qvsat,pressure
+  real :: rh,lsp,convp
+  real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
+  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
+  real :: xlon,ylat,xlonr,dzdx,dzdy
+  real :: dzdx1,dzdx2,dzdy1,dzdy2
+  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
+  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
+
+  ! Other variables:
+  ! ix,jy,kz,kzmin     loop control indices in each direction            *
+  ! kl,klp
+  ! ix1,jy1,ixp,jyp,ixm,jym 
+  integer :: ix,jy,kz,iz,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
+  !***********************************************************************
+
+  !***********************************************************************
+  ! Local constants                                                      *
+  real,parameter :: const=r_air/ga
+  !***********************************************************************
+
+  !***********************************************************************
+  ! Global variables                                                     *
+  !     from com_mod                                                     *
+  ! nxmin1,nymin1           nx-1, ny-1, respectively                     *
+  ! tt2                  2 meter temperature                             *
+  !  ps                   surface pressure                               * 
+  ! td2                  2 meter dew point                               *
+  ! height                  heights of all levels                        *
+  ! akm, bkm  coeffs.  which regulate vertical discretization of ecmwf   *
+  ! akz, bkz  model discretization coeffizients at the centre of layers  *
+  ! nuvz,nwz                vertical dimension of original ECMWF data    *
+  ! nx,ny,nz                actual dimensions of wind fields in x,y and z*
+  ! aknew,bknew model discretization coeffs. at the interpolated levels  *
+  ! uu,vv,ww [m/2]       wind components in x,y and z direction          *
+  ! qv                   specific humidity data                          *
+  ! pv (pvu)             potential vorticity                             *
+  ! rho [kg/m3]          air density                                     *
+  ! drhodz [kg/m2]       vertical air density gradient                   *
+  ! dxconst,dyconst         auxiliary variables for utransform,vtransform*
+  ! ylat0                 geographical latitude of lower left grid point *
+  ! xlon0                 geographical longitude of lower left grid point*
+  ! uupol,vvpol [m/s]   wind components in polar stereographic projection*
+  ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition     *    
+  ! lsprec [mm/h]        large scale total precipitation                 *
+  ! convprec [mm/h]      convective precipitation                        *
+  ! cloudsh                                                              *
+  !                                                                      *
+  !***********************************************************************
+
+!-----------------------------------------------------------------------------
+
+
+
+  logical :: init = .true.
+
+
+  !*************************************************************************
+  ! If verttransform is called the first time, initialize heights of the   *
+  ! z levels in meter. The heights are the heights of model levels, where  *
+  ! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
+  ! the vertical resolution in the z system is doubled. As reference point,*
+  ! the lower left corner of the grid is used.                             *
+  ! Unlike in the eta system, no difference between heights for u,v and    *
+  ! heights for w exists.                                                  *
+  !*************************************************************************
+
+
+  !      do 897 kz=1,nuvz
+  !       write (*,*) 'akz: ',akz(kz),'bkz',bkz(kz)
+  !897     continue
+
+  if (init) then
+
+  ! Search for a point with high surface pressure (i.e. not above significant topography)
+  ! Then, use this point to construct a reference z profile, to be used at all times
+  !*****************************************************************************
+
+  ! Alternatively, and to be used in GRIB2NC4, the reference point can be
+  !   provided by the user
+
+    if ( coordx_opt >= 0 .and. coordy_opt >= 0 ) then
+      print *, "Using user provided coordinates: ", coordx_opt,", ",coordy_opt
+      ixm = coordx_opt
+      jym = coordy_opt
+    else
+      print *, "Searching for coordinates"
+      do jy=0,nymin1
+        do ix=0,nxmin1
+          if (ps(ix,jy,1,n).gt.100000.) then
+            ixm=ix
+            jym=jy
+            goto 3
+          endif
+        end do
+      end do
+3     continue
+    endif
+
+
+    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
+         ps(ixm,jym,1,n))
+    pold=ps(ixm,jym,1,n)
+    height(1)=0.
+
+    do kz=2,nuvz
+      pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
+      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
+
+
+  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
+  ! upon the transformation to z levels. In order to save computer memory, this is
+  ! not done anymore in the standard version. However, this option can still be
+  ! switched on by replacing the following lines with those below, that are
+  ! currently commented out.
+  ! Note that two more changes are necessary in this subroutine below.
+  ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
+  !*****************************************************************************
+
+      if (abs(tv-tvold).gt.0.2) then
+        height(kz)= &
+             height(kz-1)+const*log(pold/pint)* &
+             (tv-tvold)/log(tv/tvold)
+      else
+        height(kz)=height(kz-1)+ &
+             const*log(pold/pint)*tv
+      endif
+
+  ! Switch on following lines to use doubled vertical resolution
+  !*************************************************************
+  !    if (abs(tv-tvold).gt.0.2) then
+  !      height((kz-1)*2)=
+  !    +      height(max((kz-2)*2,1))+const*log(pold/pint)*
+  !    +      (tv-tvold)/log(tv/tvold)
+  !    else
+  !      height((kz-1)*2)=height(max((kz-2)*2,1))+
+  !    +      const*log(pold/pint)*tv
+  !    endif
+  ! End doubled vertical resolution
+
+      tvold=tv
+      pold=pint
+    end do
+
+
+  ! Switch on following lines to use doubled vertical resolution
+  !*************************************************************
+  !  do 7 kz=3,nz-1,2
+  !    height(kz)=0.5*(height(kz-1)+height(kz+1))
+  !  height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
+  ! End doubled vertical resolution
+
+
+  ! Determine highest levels that can be within PBL
+  !************************************************
+
+    do kz=1,nz
+      if (height(kz).gt.hmixmax) then
+        nmixz=kz
+        goto 9
+      endif
+    end do
+9   continue
+
+  ! Do not repeat initialization of the Cartesian z grid
+  !*****************************************************
+
+    init=.false.
+
+  endif
+
+
+  ! Loop over the whole grid
+  !*************************
+
+  do jy=0,nymin1
+    do ix=0,nxmin1
+      tvold=tt2(ix,jy,1,n)*(1.+0.378*ew(td2(ix,jy,1,n))/ &
+           ps(ix,jy,1,n))
+      pold=ps(ix,jy,1,n)
+      uvzlev(1)=0.
+      wzlev(1)=0.
+      rhoh(1)=pold/(r_air*tvold)
+
+
+  ! Compute heights of eta levels
+  !******************************
+
+      do kz=2,nuvz
+        pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
+        tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
+        rhoh(kz)=pint/(r_air*tv)
+
+        if (abs(tv-tvold).gt.0.2) then
+          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
+               (tv-tvold)/log(tv/tvold)
+        else
+          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
+        endif
+
+        tvold=tv
+        pold=pint
+      end do
+
+
+      do kz=2,nwz-1
+        wzlev(kz)=(uvzlev(kz+1)+uvzlev(kz))/2.
+      end do
+      wzlev(nwz)=wzlev(nwz-1)+ &
+           uvzlev(nuvz)-uvzlev(nuvz-1)
+
+      uvwzlev(ix,jy,1)=0.0
+      do kz=2,nuvz
+        uvwzlev(ix,jy,kz)=uvzlev(kz)
+      end do
+
+  ! Switch on following lines to use doubled vertical resolution
+  ! Switch off the three lines above.
+  !*************************************************************
+  !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
+  !     do 23 kz=2,nwz
+  !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
+  ! End doubled vertical resolution
+
+  ! pinmconv=(h2-h1)/(p2-p1)
+
+      pinmconv(1)=(uvwzlev(ix,jy,2)-uvwzlev(ix,jy,1))/ &
+           ((aknew(2)+bknew(2)*ps(ix,jy,1,n))- &
+           (aknew(1)+bknew(1)*ps(ix,jy,1,n)))
+      do kz=2,nz-1
+        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
+             ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
+             (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
+      end do
+      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
+           ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
+           (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
+
+  ! Levels, where u,v,t and q are given
+  !************************************
+
+      uu(ix,jy,1,n)=uuh(ix,jy,1)
+      vv(ix,jy,1,n)=vvh(ix,jy,1)
+      tt(ix,jy,1,n)=tth(ix,jy,1,n)
+      qv(ix,jy,1,n)=qvh(ix,jy,1,n)
+      pv(ix,jy,1,n)=pvh(ix,jy,1)
+      rho(ix,jy,1,n)=rhoh(1)
+      uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
+      vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
+      tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
+      qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
+      pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
+      rho(ix,jy,nz,n)=rhoh(nuvz)
+      kmin=2
+      do iz=2,nz-1
+        do kz=kmin,nuvz
+          if(height(iz).gt.uvzlev(nuvz)) then
+            uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
+            vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
+            tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
+            qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
+            pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
+            rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
+            goto 30
+          endif
+          if ((height(iz).gt.uvzlev(kz-1)).and. &
+               (height(iz).le.uvzlev(kz))) then
+           dz1=height(iz)-uvzlev(kz-1)
+           dz2=uvzlev(kz)-height(iz)
+           dz=dz1+dz2
+           uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
+           vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
+           tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
+                +tth(ix,jy,kz,n)*dz1)/dz
+           qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
+                +qvh(ix,jy,kz,n)*dz1)/dz
+           pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
+           rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
+           kmin=kz
+           goto 30
+          endif
+        end do
+30      continue
+      end do
+
+
+  ! Levels, where w is given
+  !*************************
+
+      ww(ix,jy,1,n)=wwh(ix,jy,1)*pinmconv(1)
+      ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
+      kmin=2
+      do iz=2,nz
+        do kz=kmin,nwz
+          if ((height(iz).gt.wzlev(kz-1)).and. &
+               (height(iz).le.wzlev(kz))) then
+           dz1=height(iz)-wzlev(kz-1)
+           dz2=wzlev(kz)-height(iz)
+           dz=dz1+dz2
+           ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
+                +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
+           kmin=kz
+           goto 40
+          endif
+        end do
+40      continue
+      end do
+
+  ! Compute density gradients at intermediate levels
+  !*************************************************
+
+      drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
+           (height(2)-height(1))
+      do kz=2,nz-1
+        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
+             (height(kz+1)-height(kz-1))
+      end do
+      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
+
+    end do
+  end do
+
+
+  !****************************************************************
+  ! Compute slope of eta levels in windward direction and resulting
+  ! vertical wind correction
+  !****************************************************************
+
+  do jy=1,ny-2
+    do ix=1,nx-2
+
+      kmin=2
+      do iz=2,nz-1
+
+        ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
+        vi=vv(ix,jy,iz,n)*dyconst
+
+        do kz=kmin,nz
+          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
+               (height(iz).le.uvwzlev(ix,jy,kz))) then
+            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
+            dz2=uvwzlev(ix,jy,kz)-height(iz)
+            dz=dz1+dz2
+            kl=kz-1
+            klp=kz
+            kmin=kz
+            goto 47
+          endif
+        end do
+
+47      ix1=ix-1
+        jy1=jy-1
+        ixp=ix+1
+        jyp=jy+1
+
+        dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
+        dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
+        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
+
+        dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
+        dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
+        dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
+
+        ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)
+
+      end do
+
+    end do
+  end do
+
+
+  ! If north pole is in the domain, calculate wind velocities in polar
+  ! stereographic coordinates
+  !*******************************************************************
+
+  if (nglobal) then
+    do jy=int(switchnorthg)-2,nymin1
+      ylat=ylat0+real(jy)*dy
+      do ix=0,nxmin1
+        xlon=xlon0+real(ix)*dx
+        do iz=1,nz
+          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
+               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+               vvpol(ix,jy,iz,n))
+        end do
+      end do
+    end do
+
+
+    do iz=1,nz
+
+  ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+  !
+  !   AMSnauffer Nov 18 2004 Added check for case vv=0
+  !
+      xlon=xlon0+real(nx/2-1)*dx
+      xlonr=xlon*pi/180.
+      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
+           vv(nx/2-1,nymin1,iz,n)**2)
+      if (vv(nx/2-1,nymin1,iz,n).lt.0.) then
+        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
+             vv(nx/2-1,nymin1,iz,n))-xlonr
+      else if (vv(nx/2-1,nymin1,iz,n).gt.0.) then
+        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
+             vv(nx/2-1,nymin1,iz,n))-xlonr
+      else
+        ddpol=pi/2-xlonr
+      endif
+      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
+      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
+
+  ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
+      xlon=180.0
+      xlonr=xlon*pi/180.
+      ylat=90.0
+      uuaux=-ffpol*sin(xlonr+ddpol)
+      vvaux=-ffpol*cos(xlonr+ddpol)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
+
+      jy=nymin1
+      do ix=0,nxmin1
+        uupol(ix,jy,iz,n)=uupolaux
+        vvpol(ix,jy,iz,n)=vvpolaux
+      end do
+    end do
+
+
+  ! Fix: Set W at pole to the zonally averaged W of the next equator-
+  ! ward parallel of latitude
+
+  do iz=1,nz
+      wdummy=0.
+      jy=ny-2
+      do ix=0,nxmin1
+        wdummy=wdummy+ww(ix,jy,iz,n)
+      end do
+      wdummy=wdummy/real(nx)
+      jy=nymin1
+      do ix=0,nxmin1
+        ww(ix,jy,iz,n)=wdummy
+      end do
+  end do
+
+  endif
+
+
+  ! If south pole is in the domain, calculate wind velocities in polar
+  ! stereographic coordinates
+  !*******************************************************************
+
+  if (sglobal) then
+    do jy=0,int(switchsouthg)+3
+      ylat=ylat0+real(jy)*dy
+      do ix=0,nxmin1
+        xlon=xlon0+real(ix)*dx
+        do iz=1,nz
+          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
+               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+               vvpol(ix,jy,iz,n))
+        end do
+      end do
+    end do
+
+    do iz=1,nz
+
+  ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+  !
+  !   AMSnauffer Nov 18 2004 Added check for case vv=0
+  !
+      xlon=xlon0+real(nx/2-1)*dx
+      xlonr=xlon*pi/180.
+      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
+           vv(nx/2-1,0,iz,n)**2)
+      if (vv(nx/2-1,0,iz,n).lt.0.) then
+        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
+             vv(nx/2-1,0,iz,n))+xlonr
+      else if (vv(nx/2-1,0,iz,n).gt.0.) then
+        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
+             vv(nx/2-1,0,iz,n))+xlonr
+      else
+        ddpol=pi/2-xlonr
+      endif
+      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
+      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
+
+  ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
+      xlon=180.0
+      xlonr=xlon*pi/180.
+      ylat=-90.0
+      uuaux=+ffpol*sin(xlonr-ddpol)
+      vvaux=-ffpol*cos(xlonr-ddpol)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
+
+      jy=0
+      do ix=0,nxmin1
+        uupol(ix,jy,iz,n)=uupolaux
+        vvpol(ix,jy,iz,n)=vvpolaux
+      end do
+    end do
+
+
+  ! Fix: Set W at pole to the zonally averaged W of the next equator-
+  ! ward parallel of latitude
+
+    do iz=1,nz
+      wdummy=0.
+      jy=1
+      do ix=0,nxmin1
+        wdummy=wdummy+ww(ix,jy,iz,n)
+      end do
+      wdummy=wdummy/real(nx)
+      jy=0
+      do ix=0,nxmin1
+        ww(ix,jy,iz,n)=wdummy
+      end do
+    end do
+  endif
+
+
+  !write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
+  !   create a cloud and rainout/washout field, clouds occur where rh>80%
+  !   total cloudheight is stored at level 0
+  do jy=0,nymin1
+    do ix=0,nxmin1
+      rain_cloud_above=0
+      lsp=lsprec(ix,jy,1,n)
+      convp=convprec(ix,jy,1,n)
+      cloudsh(ix,jy,n)=0
+      do kz_inv=1,nz-1
+         kz=nz-kz_inv+1
+         pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+         rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+         clouds(ix,jy,kz,n)=0
+         if (rh.gt.0.8) then ! in cloud
+            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+               rain_cloud_above=1
+               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
+                    height(kz)-height(kz-1)
+               if (lsp.ge.convp) then
+                  clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
+               else
+                  clouds(ix,jy,kz,n)=2 ! convp dominated rainout
+               endif
+            else ! no precipitation
+                  clouds(ix,jy,kz,n)=1 ! cloud
+            endif
+         else ! no cloud
+            if (rain_cloud_above.eq.1) then ! scavenging
+               if (lsp.ge.convp) then
+                  clouds(ix,jy,kz,n)=5 ! lsp dominated washout
+               else
+                  clouds(ix,jy,kz,n)=4 ! convp dominated washout
+               endif
+            endif
+         endif
+      end do
+    end do
+  end do
+
+  !do 102 kz=1,nuvz
+  !write(an,'(i02)') kz+10
+  !write(*,*) nuvz,nymin1,nxmin1,'--',an,'--'
+  !open(4,file='/nilu_wrk2/sec/cloudtest/cloud'//an,form='formatted')
+  !do 101 jy=0,nymin1
+  !    write(4,*) (clouds(ix,jy,kz,n),ix=1,nxmin1)
+  !101   continue
+  ! close(4)
+  !102   continue
+
+  ! open(4,file='/nilu_wrk2/sec/cloudtest/height',form='formatted')
+  ! do 103 jy=0,nymin1
+  !     write (4,*)
+  !+       (height(kz),kz=1,nuvz)
+  !103   continue
+  ! close(4)
+
+  !open(4,file='/nilu_wrk2/sec/cloudtest/p',form='formatted')
+  ! do 104 jy=0,nymin1
+  !     write (4,*)
+  !+       (r_air*tt(ix,jy,1,n)*rho(ix,jy,1,n),ix=1,nxmin1)
+  !104   continue
+  ! close(4)
+end subroutine verttransform_grib2nc4_ecmwf
Index: flexpart_code/grib2nc4/verttransform_grib2nc4_gfs.f90
===================================================================
--- flexpart_code/grib2nc4/verttransform_grib2nc4_gfs.f90	(revision 8624a75e20c0cc6dff42b4eb97182a138f9aa1e7)
+++ flexpart_code/grib2nc4/verttransform_grib2nc4_gfs.f90	(revision 8624a75e20c0cc6dff42b4eb97182a138f9aa1e7)
@@ -0,0 +1,680 @@
+!**********************************************************************
+! Copyright 1998,1999,2000,2001,2002,2005,2007,2008,2009,2010         *
+! Andreas Stohl, Petra Seibert, A. Frank, Gerhard Wotawa,             *
+! Caroline Forster, Sabine Eckhardt, John Burkhart, Harald Sodemann   *
+!                                                                     *
+! This file is part of FLEXPART.                                      *
+!                                                                     *
+! FLEXPART is free software: you can redistribute it and/or modify    *
+! it under the terms of the GNU General Public License as published by*
+! the Free Software Foundation, either version 3 of the License, or   *
+! (at your option) any later version.                                 *
+!                                                                     *
+! FLEXPART is distributed in the hope that it will be useful,         *
+! but WITHOUT ANY WARRANTY; without even the implied warranty of      *
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the       *
+! GNU General Public License for more details.                        *
+!                                                                     *
+! You should have received a copy of the GNU General Public License   *
+! along with FLEXPART.  If not, see <http://www.gnu.org/licenses/>.   *
+!**********************************************************************
+
+subroutine verttransform_grib2nc4_gfs(n,uuh,vvh,wwh,pvh,coordx_opt,coordy_opt)
+  !                         i  i   i   i   i
+  !*****************************************************************************
+  !                                                                            *
+  !     This subroutine transforms temperature, dew point temperature and      *
+  !     wind components from eta to meter coordinates.                         *
+  !     The vertical wind component is transformed from Pa/s to m/s using      *
+  !     the conversion factor pinmconv.                                        *
+  !     In addition, this routine calculates vertical density gradients        *
+  !     needed for the parameterization of the turbulent velocities.           *
+  !                                                                            *
+  !     Author: A. Stohl, G. Wotawa                                            *
+  !                                                                            *
+  !     12 August 1996                                                         *
+  !     Update: 16 January 1998                                                *
+  !                                                                            *
+  !     Major update: 17 February 1999                                         *
+  !     by G. Wotawa                                                           *
+  !     CHANGE 17/11/2005 Caroline Forster, NCEP GFS version                   *
+  !                                                                            *
+  !   - Vertical levels for u, v and w are put together                        *
+  !   - Slope correction for vertical velocity: Modification of calculation    *
+  !     procedure                                                              *
+  !                                                                            *
+  !*****************************************************************************
+  !  Changes, Bernd C. Krueger, Feb. 2001:
+  !   Variables tth and qvh (on eta coordinates) from common block
+  !*****************************************************************************
+  !  Changes Arnold, D. and Morton, D. (2015):                                 *
+  !   -- description of local and common variables                             *
+  !***************************************************************************** 
+  !  July 2017 - M. Harustak                                                   *
+  !  In order to be used by GRIB2NC4, this routine has been modified so that   *
+  !  it has the possibility to get two additional arguments, coordx_opt and    *
+  !  coordy_opt. These two are used to define the location where the FLEXPART  *
+  !  model levels are calculated and used in GRIB2NC4 to obtaine the met data  *
+  !                                                                            *
+  !  The routine works normally if -1 (no value) is passed                     *
+  !*****************************************************************************
+
+
+  use par_mod
+  use com_mod
+  use cmapf_mod
+
+  implicit none
+  !***********************************************************************
+  ! Subroutine Parameters:                                               *
+  !    input                                                             *
+  ! n                   temporal index for meteorological fields (1 to 3)*  
+  ! uuh,vvh, wwh        wind components in ECMWF model levels            *  
+  ! pvh                 potential vorticity in ECMWF model levels        *
+  integer :: n
+  real :: uuh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real :: vvh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real :: pvh(0:nxmax-1,0:nymax-1,nuvzmax)
+  real :: wwh(0:nxmax-1,0:nymax-1,nwzmax)
+  integer :: coordx_opt, coordy_opt
+  !***********************************************************************
+
+  !***********************************************************************
+  ! Local variables                                                      *
+  !                                                                      *
+  !  ew                   subroutine/function to calculate saturation    *  
+  !                         water vaport for a given air temperature     *
+  ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition     *
+  ! rain_cloud_above [0,1]  whether there is a raining cloud or not      *  
+  ! kz_inv  inverted indez for kz                                        *
+  ! f_qvsat  Saturation water vapor specific humidity (kg/kg)            *
+  ! pressure                                                             *
+  ! rh       relative humidity                                           *    
+  ! lsp      large scale precip at one grid cell                         *
+  ! convp    convective precip at one grid cell                          *
+  ! uvzlev    height of the eta half-levels                              *
+  ! rhoh     density in model levels                                     *
+  ! pinmconv conversion factor                                           *
+  ! pint     pressure on model levels (using akz, bkz, ps)               *
+  ! tv       virtual temperature                                         *
+  ! tvold,pold      dummy variables to keep previous value               *
+  ! dz1, dz2, dz  differences heights model levels, pressure levels      *    
+  ! ui, vi                                                               *   
+  ! xlon,xlat                                                            *
+  ! xlonr         xlon*pi/180.                                           * 
+  ! dzdx,dzdy     slope corrections                                      *
+  ! dzdx1,dzdx2,dzdy1   slope corrections in each direction              *
+  ! uuaux,vvaux,uupolaux,vvpolaux  auxiliary variables for polar sterero.*
+  ! ddpol,ffpol  for special treatment of poles                          *
+  ! wdummy                                                               *
+  ! wzlev                                                                *
+  ! uvwzlev                                                              *
+
+  integer :: rain_cloud_above,kz_inv
+  real :: f_qvsat,pressure
+  real :: rh,lsp,convp
+  real :: uvzlev(nuvzmax),rhoh(nuvzmax),pinmconv(nzmax)
+  real :: ew,pint,tv,tvold,pold,dz1,dz2,dz,ui,vi
+  real :: xlon,ylat,xlonr,dzdx,dzdy
+  real :: dzdx1,dzdx2,dzdy1,dzdy2
+  real :: uuaux,vvaux,uupolaux,vvpolaux,ddpol,ffpol,wdummy
+  real :: wzlev(nwzmax),uvwzlev(0:nxmax-1,0:nymax-1,nzmax)
+
+  ! NCEP version
+  integer :: llev, i
+
+  ! Other variables:
+  ! ix,jy,kz,kzmin     loop control indices in each direction            *
+  ! kl,klp
+  ! ix1,jy1,ixp,jyp,ixm,jym 
+  integer :: ix,jy,kz,iz,kmin,kl,klp,ix1,jy1,ixp,jyp,ixm,jym
+  !***********************************************************************
+
+  !***********************************************************************
+  ! Local constants                                                      *
+  real,parameter :: const=r_air/ga
+  !***********************************************************************
+
+  !***********************************************************************
+  ! Global variables                                                     *
+  !     from par_mod and com_mod                                         *
+  ! nxmin1,nymin1           nx-1, ny-1, respectively                     *
+  ! tt2                  2 meter temperature                             *
+  !  ps                   surface pressure                               * 
+  ! td2                  2 meter dew point                               *
+  ! height                  heights of all levels                        *
+  ! akm, bkm  coeffs.  which regulate vertical discretization of ecmwf   *
+  ! akz, bkz  model discretization coeffizients at the centre of layers  *
+  ! tv    virtual temperature                                            *
+  ! nuvz,nwz                vertical dimension of original ECMWF data    *
+  ! nx,ny,nz                actual dimensions of wind fields in x,y and z*
+  ! aknew,bknew model discretization coeffs. at the interpolated levels  *
+  ! uu,vv,ww [m/2]       wind components in x,y and z direction          *
+  ! qv                   specific humidity data                          *
+  ! pv (pvu)             potential vorticity                             *
+  ! rho [kg/m3]          air density                                     *
+  ! drhodz [kg/m2]       vertical air density gradient                   *
+  ! dxconst,dyconst         auxiliary variables for utransform,vtransform*
+  ! ylat0                 geographical latitude of lower left grid point *
+  ! xlon0                 geographical longitude of lower left grid point*
+  ! uupol,vvpol [m/s]   wind components in polar stereographic projection*
+  ! clouds(0:nxmax,0:nymax,0:nzmax,2) cloud field for wet deposition     *    
+  ! lsprec [mm/h]        large scale total precipitation                 *
+  ! convprec [mm/h]      convective precipitation                        *
+  ! cloudsh                                                              *
+  !                                                                      *
+  !***********************************************************************
+
+!-----------------------------------------------------------------------------
+
+
+  logical :: init = .true.
+
+
+  !*************************************************************************
+  ! If verttransform is called the first time, initialize heights of the   *
+  ! z levels in meter. The heights are the heights of model levels, where  *
+  ! u,v,T and qv are given, and of the interfaces, where w is given. So,   *
+  ! the vertical resolution in the z system is doubled. As reference point,*
+  ! the lower left corner of the grid is used.                             *
+  ! Unlike in the eta system, no difference between heights for u,v and    *
+  ! heights for w exists.                                                  *
+  !*************************************************************************
+
+  if (init) then
+
+  ! Search for a point with high surface pressure (i.e. not above significant topography)
+  ! Then, use this point to construct a reference z profile, to be used at all times
+  !*****************************************************************************
+
+
+    if ( coordx_opt >= 0 .and. coordy_opt>=0 ) then
+      print *, "Using user provided coordinates: ", coordx_opt,", ",coordy_opt
+      ixm = coordx_opt
+      jym = coordy_opt
+    else
+      do jy=0,nymin1
+        do ix=0,nxmin1
+          if (ps(ix,jy,1,n).gt.100000.) then
+            ixm=ix
+            jym=jy
+            goto 3
+          endif
+        end do
+      end do
+3     continue
+    endif
+
+
+    tvold=tt2(ixm,jym,1,n)*(1.+0.378*ew(td2(ixm,jym,1,n))/ &
+         ps(ixm,jym,1,n))
+    pold=ps(ixm,jym,1,n)
+    height(1)=0.
+
+    do kz=2,nuvz
+      pint=akz(kz)+bkz(kz)*ps(ixm,jym,1,n)
+      tv=tth(ixm,jym,kz,n)*(1.+0.608*qvh(ixm,jym,kz,n))
+
+
+  ! NOTE: In FLEXPART versions up to 4.0, the number of model levels was doubled
+  ! upon the transformation to z levels. In order to save computer memory, this is
+  ! not done anymore in the standard version. However, this option can still be
+  ! switched on by replacing the following lines with those below, that are
+  ! currently commented out.
+  ! Note that two more changes are necessary in this subroutine below.
+  ! One change is also necessary in gridcheck.f, and another one in verttransform_nests.
+  !*****************************************************************************
+
+      if (abs(tv-tvold).gt.0.2) then
+        height(kz)= &
+             height(kz-1)+const*log(pold/pint)* &
+             (tv-tvold)/log(tv/tvold)
+      else
+        height(kz)=height(kz-1)+ &
+             const*log(pold/pint)*tv
+      endif
+
+  ! Switch on following lines to use doubled vertical resolution
+  !*************************************************************
+  !    if (abs(tv-tvold).gt.0.2) then
+  !      height((kz-1)*2)=
+  !    +      height(max((kz-2)*2,1))+const*log(pold/pint)*
+  !    +      (tv-tvold)/log(tv/tvold)
+  !    else
+  !      height((kz-1)*2)=height(max((kz-2)*2,1))+
+  !    +      const*log(pold/pint)*tv
+  !    endif
+  ! End doubled vertical resolution
+
+      tvold=tv
+      pold=pint
+    end do
+
+
+  ! Switch on following lines to use doubled vertical resolution
+  !*************************************************************
+  !  do 7 kz=3,nz-1,2
+  !    height(kz)=0.5*(height(kz-1)+height(kz+1))
+  !  height(nz)=height(nz-1)+height(nz-1)-height(nz-2)
+  ! End doubled vertical resolution
+
+
+  ! Determine highest levels that can be within PBL
+  !************************************************
+
+    do kz=1,nz
+      if (height(kz).gt.hmixmax) then
+        nmixz=kz
+        goto 9
+      endif
+    end do
+9   continue
+
+  ! Do not repeat initialization of the Cartesian z grid
+  !*****************************************************
+
+    init=.false.
+
+  endif
+
+
+  ! Loop over the whole grid
+  !*************************
+
+  do jy=0,nymin1
+    do ix=0,nxmin1
+
+  ! NCEP version: find first level above ground
+      llev = 0
+      do i=1,nuvz
+       if (ps(ix,jy,1,n).lt.akz(i)) llev=i
+      end do
+       llev = llev+1
+       if (llev.gt.nuvz-2) llev = nuvz-2
+  !     if (llev.eq.nuvz-2) write(*,*) 'verttransform
+  !    +WARNING: LLEV eq NUZV-2'
+  ! NCEP version
+
+
+  ! compute height of pressure levels above ground
+  !***********************************************
+
+      tvold=tth(ix,jy,llev,n)*(1.+0.608*qvh(ix,jy,llev,n))
+      pold=akz(llev)
+      uvzlev(llev)=0.
+      wzlev(llev)=0.
+      uvwzlev(ix,jy,llev)=0.
+      rhoh(llev)=pold/(r_air*tvold)
+
+      do kz=llev+1,nuvz
+        pint=akz(kz)+bkz(kz)*ps(ix,jy,1,n)
+        tv=tth(ix,jy,kz,n)*(1.+0.608*qvh(ix,jy,kz,n))
+        rhoh(kz)=pint/(r_air*tv)
+
+        if (abs(tv-tvold).gt.0.2) then
+          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)* &
+               (tv-tvold)/log(tv/tvold)
+        else
+          uvzlev(kz)=uvzlev(kz-1)+const*log(pold/pint)*tv
+        endif
+        wzlev(kz)=uvzlev(kz)
+        uvwzlev(ix,jy,kz)=uvzlev(kz)
+
+        tvold=tv
+        pold=pint
+      end do
+
+
+  ! Switch on following lines to use doubled vertical resolution
+  ! Switch off the three lines above.
+  !*************************************************************
+  !22          uvwzlev(ix,jy,(kz-1)*2)=uvzlev(kz)
+  !     do 23 kz=2,nwz
+  !23          uvwzlev(ix,jy,(kz-1)*2+1)=wzlev(kz)
+  ! End doubled vertical resolution
+
+  ! pinmconv=(h2-h1)/(p2-p1)
+
+      pinmconv(llev)=(uvwzlev(ix,jy,llev+1)-uvwzlev(ix,jy,llev))/ &
+           ((aknew(llev+1)+bknew(llev+1)*ps(ix,jy,1,n))- &
+           (aknew(llev)+bknew(llev)*ps(ix,jy,1,n)))
+      do kz=llev+1,nz-1
+        pinmconv(kz)=(uvwzlev(ix,jy,kz+1)-uvwzlev(ix,jy,kz-1))/ &
+             ((aknew(kz+1)+bknew(kz+1)*ps(ix,jy,1,n))- &
+             (aknew(kz-1)+bknew(kz-1)*ps(ix,jy,1,n)))
+      end do
+      pinmconv(nz)=(uvwzlev(ix,jy,nz)-uvwzlev(ix,jy,nz-1))/ &
+           ((aknew(nz)+bknew(nz)*ps(ix,jy,1,n))- &
+           (aknew(nz-1)+bknew(nz-1)*ps(ix,jy,1,n)))
+
+
+  ! Levels, where u,v,t and q are given
+  !************************************
+
+      uu(ix,jy,1,n)=uuh(ix,jy,llev)
+      vv(ix,jy,1,n)=vvh(ix,jy,llev)
+      tt(ix,jy,1,n)=tth(ix,jy,llev,n)
+      qv(ix,jy,1,n)=qvh(ix,jy,llev,n)
+      pv(ix,jy,1,n)=pvh(ix,jy,llev)
+      rho(ix,jy,1,n)=rhoh(llev)
+      pplev(ix,jy,1,n)=akz(llev)
+      uu(ix,jy,nz,n)=uuh(ix,jy,nuvz)
+      vv(ix,jy,nz,n)=vvh(ix,jy,nuvz)
+      tt(ix,jy,nz,n)=tth(ix,jy,nuvz,n)
+      qv(ix,jy,nz,n)=qvh(ix,jy,nuvz,n)
+      pv(ix,jy,nz,n)=pvh(ix,jy,nuvz)
+      rho(ix,jy,nz,n)=rhoh(nuvz)
+      pplev(ix,jy,nz,n)=akz(nuvz)
+      kmin=llev+1
+      do iz=2,nz-1
+        do kz=kmin,nuvz
+          if(height(iz).gt.uvzlev(nuvz)) then
+            uu(ix,jy,iz,n)=uu(ix,jy,nz,n)
+            vv(ix,jy,iz,n)=vv(ix,jy,nz,n)
+            tt(ix,jy,iz,n)=tt(ix,jy,nz,n)
+            qv(ix,jy,iz,n)=qv(ix,jy,nz,n)
+            pv(ix,jy,iz,n)=pv(ix,jy,nz,n)
+            rho(ix,jy,iz,n)=rho(ix,jy,nz,n)
+            pplev(ix,jy,iz,n)=pplev(ix,jy,nz,n)
+            goto 30
+          endif
+          if ((height(iz).gt.uvzlev(kz-1)).and. &
+               (height(iz).le.uvzlev(kz))) then
+           dz1=height(iz)-uvzlev(kz-1)
+           dz2=uvzlev(kz)-height(iz)
+           dz=dz1+dz2
+           uu(ix,jy,iz,n)=(uuh(ix,jy,kz-1)*dz2+uuh(ix,jy,kz)*dz1)/dz
+           vv(ix,jy,iz,n)=(vvh(ix,jy,kz-1)*dz2+vvh(ix,jy,kz)*dz1)/dz
+           tt(ix,jy,iz,n)=(tth(ix,jy,kz-1,n)*dz2 &
+                +tth(ix,jy,kz,n)*dz1)/dz
+           qv(ix,jy,iz,n)=(qvh(ix,jy,kz-1,n)*dz2 &
+                +qvh(ix,jy,kz,n)*dz1)/dz
+           pv(ix,jy,iz,n)=(pvh(ix,jy,kz-1)*dz2+pvh(ix,jy,kz)*dz1)/dz
+           rho(ix,jy,iz,n)=(rhoh(kz-1)*dz2+rhoh(kz)*dz1)/dz
+           pplev(ix,jy,iz,n)=(akz(kz-1)*dz2+akz(kz)*dz1)/dz
+          endif
+        end do
+30      continue
+      end do
+
+
+  ! Levels, where w is given
+  !*************************
+
+      ww(ix,jy,1,n)=wwh(ix,jy,llev)*pinmconv(llev)
+      ww(ix,jy,nz,n)=wwh(ix,jy,nwz)*pinmconv(nz)
+      kmin=llev+1
+      do iz=2,nz
+        do kz=kmin,nwz
+          if ((height(iz).gt.wzlev(kz-1)).and. &
+               (height(iz).le.wzlev(kz))) then
+           dz1=height(iz)-wzlev(kz-1)
+           dz2=wzlev(kz)-height(iz)
+           dz=dz1+dz2
+           ww(ix,jy,iz,n)=(wwh(ix,jy,kz-1)*pinmconv(kz-1)*dz2 &
+                +wwh(ix,jy,kz)*pinmconv(kz)*dz1)/dz
+
+          endif
+        end do
+      end do
+
+
+  ! Compute density gradients at intermediate levels
+  !*************************************************
+
+      drhodz(ix,jy,1,n)=(rho(ix,jy,2,n)-rho(ix,jy,1,n))/ &
+           (height(2)-height(1))
+      do kz=2,nz-1
+        drhodz(ix,jy,kz,n)=(rho(ix,jy,kz+1,n)-rho(ix,jy,kz-1,n))/ &
+             (height(kz+1)-height(kz-1))
+      end do
+      drhodz(ix,jy,nz,n)=drhodz(ix,jy,nz-1,n)
+
+    end do
+  end do
+
+
+  !****************************************************************
+  ! Compute slope of eta levels in windward direction and resulting
+  ! vertical wind correction
+  !****************************************************************
+
+  do jy=1,ny-2
+    do ix=1,nx-2
+
+  ! NCEP version: find first level above ground
+      llev = 0
+      do i=1,nuvz
+       if (ps(ix,jy,1,n).lt.akz(i)) llev=i
+      end do
+       llev = llev+1
+       if (llev.gt.nuvz-2) llev = nuvz-2
+  !     if (llev.eq.nuvz-2) write(*,*) 'verttransform
+  !    +WARNING: LLEV eq NUZV-2'
+  ! NCEP version
+
+      kmin=llev+1
+      do iz=2,nz-1
+
+        ui=uu(ix,jy,iz,n)*dxconst/cos((real(jy)*dy+ylat0)*pi180)
+        vi=vv(ix,jy,iz,n)*dyconst
+
+        do kz=kmin,nz
+          if ((height(iz).gt.uvwzlev(ix,jy,kz-1)).and. &
+               (height(iz).le.uvwzlev(ix,jy,kz))) then
+            dz1=height(iz)-uvwzlev(ix,jy,kz-1)
+            dz2=uvwzlev(ix,jy,kz)-height(iz)
+            dz=dz1+dz2
+            kl=kz-1
+            klp=kz
+            goto 47
+          endif
+        end do
+
+47      ix1=ix-1
+        jy1=jy-1
+        ixp=ix+1
+        jyp=jy+1
+
+        dzdx1=(uvwzlev(ixp,jy,kl)-uvwzlev(ix1,jy,kl))/2.
+        dzdx2=(uvwzlev(ixp,jy,klp)-uvwzlev(ix1,jy,klp))/2.
+        dzdx=(dzdx1*dz2+dzdx2*dz1)/dz
+
+        dzdy1=(uvwzlev(ix,jyp,kl)-uvwzlev(ix,jy1,kl))/2.
+        dzdy2=(uvwzlev(ix,jyp,klp)-uvwzlev(ix,jy1,klp))/2.
+        dzdy=(dzdy1*dz2+dzdy2*dz1)/dz
+
+        ww(ix,jy,iz,n)=ww(ix,jy,iz,n)+(dzdx*ui+dzdy*vi)
+
+      end do
+
+    end do
+  end do
+
+
+  ! If north pole is in the domain, calculate wind velocities in polar
+  ! stereographic coordinates
+  !*******************************************************************
+
+  if (nglobal) then
+    do jy=int(switchnorthg)-2,nymin1
+      ylat=ylat0+real(jy)*dy
+      do ix=0,nxmin1
+        xlon=xlon0+real(ix)*dx
+        do iz=1,nz
+          call cc2gll(northpolemap,ylat,xlon,uu(ix,jy,iz,n), &
+               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+               vvpol(ix,jy,iz,n))
+        end do
+      end do
+    end do
+
+
+    do iz=1,nz
+
+  ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+      xlon=xlon0+real(nx/2-1)*dx
+      xlonr=xlon*pi/180.
+      ffpol=sqrt(uu(nx/2-1,nymin1,iz,n)**2+ &
+           vv(nx/2-1,nymin1,iz,n)**2)
+      if(vv(nx/2-1,nymin1,iz,n).lt.0.) then
+        ddpol=atan(uu(nx/2-1,nymin1,iz,n)/ &
+             vv(nx/2-1,nymin1,iz,n))-xlonr
+      elseif (vv(nx/2-1,nymin1,iz,n).gt.0.) then
+        ddpol=pi+atan(uu(nx/2-1,nymin1,iz,n)/ &
+             vv(nx/2-1,nymin1,iz,n))-xlonr
+      else
+        ddpol=pi/2-xlonr
+      endif
+      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
+      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
+
+  ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
+      xlon=180.0
+      xlonr=xlon*pi/180.
+      ylat=90.0
+      uuaux=-ffpol*sin(xlonr+ddpol)
+      vvaux=-ffpol*cos(xlonr+ddpol)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
+
+      jy=nymin1
+      do ix=0,nxmin1
+        uupol(ix,jy,iz,n)=uupolaux
+        vvpol(ix,jy,iz,n)=vvpolaux
+      end do
+    end do
+
+
+  ! Fix: Set W at pole to the zonally averaged W of the next equator-
+  ! ward parallel of latitude
+
+  do iz=1,nz
+      wdummy=0.
+      jy=ny-2
+      do ix=0,nxmin1
+        wdummy=wdummy+ww(ix,jy,iz,n)
+      end do
+      wdummy=wdummy/real(nx)
+      jy=nymin1
+      do ix=0,nxmin1
+        ww(ix,jy,iz,n)=wdummy
+      end do
+  end do
+
+  endif
+
+
+  ! If south pole is in the domain, calculate wind velocities in polar
+  ! stereographic coordinates
+  !*******************************************************************
+
+  if (sglobal) then
+    do jy=0,int(switchsouthg)+3
+      ylat=ylat0+real(jy)*dy
+      do ix=0,nxmin1
+        xlon=xlon0+real(ix)*dx
+        do iz=1,nz
+          call cc2gll(southpolemap,ylat,xlon,uu(ix,jy,iz,n), &
+               vv(ix,jy,iz,n),uupol(ix,jy,iz,n), &
+               vvpol(ix,jy,iz,n))
+        end do
+      end do
+    end do
+
+    do iz=1,nz
+
+  ! CALCULATE FFPOL, DDPOL FOR CENTRAL GRID POINT
+      xlon=xlon0+real(nx/2-1)*dx
+      xlonr=xlon*pi/180.
+      ffpol=sqrt(uu(nx/2-1,0,iz,n)**2+ &
+           vv(nx/2-1,0,iz,n)**2)
+      if(vv(nx/2-1,0,iz,n).lt.0.) then
+        ddpol=atan(uu(nx/2-1,0,iz,n)/ &
+             vv(nx/2-1,0,iz,n))+xlonr
+      elseif (vv(nx/2-1,0,iz,n).gt.0.) then
+        ddpol=pi+atan(uu(nx/2-1,0,iz,n)/ &
+             vv(nx/2-1,0,iz,n))-xlonr
+      else
+        ddpol=pi/2-xlonr
+      endif
+      if(ddpol.lt.0.) ddpol=2.0*pi+ddpol
+      if(ddpol.gt.2.0*pi) ddpol=ddpol-2.0*pi
+
+  ! CALCULATE U,V FOR 180 DEG, TRANSFORM TO POLAR STEREOGRAPHIC GRID
+      xlon=180.0
+      xlonr=xlon*pi/180.
+      ylat=-90.0
+      uuaux=+ffpol*sin(xlonr-ddpol)
+      vvaux=-ffpol*cos(xlonr-ddpol)
+      call cc2gll(northpolemap,ylat,xlon,uuaux,vvaux,uupolaux, &
+           vvpolaux)
+
+      jy=0
+      do ix=0,nxmin1
+        uupol(ix,jy,iz,n)=uupolaux
+        vvpol(ix,jy,iz,n)=vvpolaux
+      end do
+    end do
+
+
+  ! Fix: Set W at pole to the zonally averaged W of the next equator-
+  ! ward parallel of latitude
+
+    do iz=1,nz
+      wdummy=0.
+      jy=1
+      do ix=0,nxmin1
+        wdummy=wdummy+ww(ix,jy,iz,n)
+      end do
+      wdummy=wdummy/real(nx)
+      jy=0
+      do ix=0,nxmin1
+        ww(ix,jy,iz,n)=wdummy
+      end do
+    end do
+  endif
+
+
+  !   write (*,*) 'initializing clouds, n:',n,nymin1,nxmin1,nz
+  !   create a cloud and rainout/washout field, clouds occur where rh>80%
+  !   total cloudheight is stored at level 0
+  do jy=0,nymin1
+    do ix=0,nxmin1
+      rain_cloud_above=0
+      lsp=lsprec(ix,jy,1,n)
+      convp=convprec(ix,jy,1,n)
+      cloudsh(ix,jy,n)=0
+      do kz_inv=1,nz-1
+         kz=nz-kz_inv+1
+         pressure=rho(ix,jy,kz,n)*r_air*tt(ix,jy,kz,n)
+         rh=qv(ix,jy,kz,n)/f_qvsat(pressure,tt(ix,jy,kz,n))
+         clouds(ix,jy,kz,n)=0
+         if (rh.gt.0.8) then ! in cloud
+            if ((lsp.gt.0.01).or.(convp.gt.0.01)) then ! cloud and precipitation
+               rain_cloud_above=1
+               cloudsh(ix,jy,n)=cloudsh(ix,jy,n)+ &
+                    height(kz)-height(kz-1)
+               if (lsp.ge.convp) then
+                  clouds(ix,jy,kz,n)=3 ! lsp dominated rainout
+               else
+                  clouds(ix,jy,kz,n)=2 ! convp dominated rainout
+               endif
+            else ! no precipitation
+                  clouds(ix,jy,kz,n)=1 ! cloud
+            endif
+         else ! no cloud
+            if (rain_cloud_above.eq.1) then ! scavenging
+               if (lsp.ge.convp) then
+                  clouds(ix,jy,kz,n)=5 ! lsp dominated washout
+               else
+                  clouds(ix,jy,kz,n)=4 ! convp dominated washout
+               endif
+            endif
+         endif
+      end do
+    end do
+  end do
+
+
+end subroutine verttransform_grib2nc4_gfs
