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)
