anchor.diff:      integer :: dimtwo(2)
       integer :: dimtwo_time(3)
       !> variable IDs
  -    integer :: lat_id,lon_id,time_id
  +    integer :: lat_id,lon_id,time_id, lev_id
       !> array to store ALL var ids
       integer :: varid(nvars)
   
  
positive.code:             varid(13)=temp_id
  
              call check( nf90_def_var(ncid, "z", NF90_REAL, dimids(1:3), temp_id), trim(err)//"z" )
              call check( nf90_put_att(ncid,temp_id,"standard_name","height"))
              call check( nf90_put_att(ncid,temp_id,"long_name","Model level height (ASL)"))
              call check( nf90_put_att(ncid,temp_id,"units","m"))
              call check( nf90_put_att(ncid,temp_id,"axis","Z"))
              ! call check( nf90_put_att(ncid,temp_id,"coordinates","lon lat"))
              varid(20)=temp_id
  
              call check( nf90_def_var(ncid, "rho", NF90_REAL, dimids, temp_id), trim(err)//"rho" )
  
negative.code:         call check( nf90_put_att(ncid,time_id,"standard_name","time"))
          call check( nf90_put_att(ncid,time_id,"UTCoffset","0"))
  
          call check( nf90_def_var(ncid, "lev", NF90_REAL, dimids(3), lev_id), trim(err)//"lev" )
          call check( nf90_put_att(ncid,lev_id,"standard_name","height"))
          call check( nf90_put_att(ncid,lev_id,"long_name","height above ground level"))
          call check( nf90_put_att(ncid,lev_id,"description","model level height (AGL)"))
  
          if (calendar==GREGORIAN) then
              call check( nf90_put_att(ncid,time_id,"long_name","modified Julian Day"))
  